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This review explains the relationship between density functional theory and strongly 
correlated models using the simplest possible example, the two-site Hubbard model. 

The relationship to traditional quantum chemistry is included. Even in this elemen¬ 
tary example, where the exact ground-state energy and site occupations can be found 
analytically, there is much to be explained in terms of the underlying logic and aims 
of Density Functional Theory. Although the usual solution is analytic, the density 
functional is given only implicitly. We overcome this difficulty using the Levy-Lieb 
construction to create a parametrization of the exact function with negligible errors. The 
symmetric case is most commonly studied, but we find a rich variation in behavior by 
including asymmetry, as strong correlation physics vies with charge-transfer effects. We 
explore the behavior of the gap and the many-body Green’s function, demonstrating 
the ‘failure’ of the Kohn-Sham method to reproduce the fundamental gap. We perform 
benchmark calculations of the occupation and components of the KS potentials, the 
correlation kinetic energies, and the adiabatic connection. We test several approximate 
functionals (restricted and unrestricted Hartree-Fock and Bethe Ansatz Local Density 
Approximation) to show their successes and limitations. We also discuss and illustrate the 
concept of the derivative discontinuity. Useful appendices include analytic expressions for 
Density Functional energy components, several limits of the exact functional (weak- and 
strong-coupling, symmetric and asymmetric), the Kohn-Sham hopping energy functional 
for 3 sites, various adiabatic connection results, proofs of exact conditions for this model, 
and the origin of the Hubbard model from a minimal basis model for stretched H 2 . 


PACS numbers: 71.15.Mb, 71.10.Fd, 71.27.+a 


1. INTRODUCTION 

In condensed matter, the world of electronic structure 
theory can be divided into two camps: the weakly and the 
strongly correlated. Weakly correlated solids are almost 
always treated with density-functional methods as a start¬ 
ing point for ground-state properties (Burke, 2012; Burke 
and Wagner, 2013; Capelle, 2006; Dreizler and Gross, 
1990; Kohn, 1999). Many-body (MB) approximations 
such as GW might then be applied to find properties of 
the quasi-particle spectrum, such as the gap(Aryasetiawan 
and Gunnarsson, 1998; Pollehn et a/., 1998; Verdozzi et a/., 
1995). This approach is ‘first-principles’, in the sense that 
it uses the real-space Hamiltonian for the electrons in the 
field of the nuclei, and produces a converged result that is 
independent of the basis set, once a sufficiently large basis 
set is used. Density functional theory (DFT) is known 
to be exact in principle, but the usual approximations 
often fail when correlations become strong(Cohen et a/., 
2008b). 

On the other hand, strongly correlated systems are most 


often treated via lattice Hamiltonians with relatively few 
parameters (Dagotto, 1994; Korepin and Essler, 1994). 
These simplified Hamiltonians can be easier to deal with, 
especially when correlations are strong(Dagotto, 1994; 
Essler et al.^ 1992). Even approximate solutions to such 
Hamiltonians can yield insight into the physics, especially 
for extended systems (Solovyev, 2008). However, such 
Hamiltonians can rarely be unambiguously derived from 
a first-principles starting point, making it difficult (if 
not impossible) to say how accurate such solutions are 
quantitatively or to improve on that accuracy. Moreover, 
methods that yield approximate Green’s functions are 
often more focused on response properties or thermal 
properties rather than on total energies in the ground- 
state. 

On the other hand, the ground-state energy of electrons 
plays a much more crucial role in chemical and material 
science applications (Martin, 2004; Parr and Yang, 1989). 
Very small energy differences determine geometries and 
sometimes qualitative properties, such as the nature of a 
transition state in a chemical reaction(Feller and Peterson, 
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2007; Helgaker et al, 2004; Lee and Scuseria, 1995) or 
where a molecule is adsorbed on a surface (B. Hammer 
and N0rskov, 1996; Over et a/., 2000). An error of 0.05 
eV changes a reaction rate by a factor of 5 at room 
temperature. Thus quantum chemical development has 
focused on extracting extremely accurate energies for 
the ground and other eigenstates(Friesner, 2005; Head- 
Gordon, 1996; Kohn et al, 1996; Schaefer, 2012; Yang 
et al.^ 2014). This is routinely achieved for molecules 
using coupled-cluster methods (CCSD(T)) and reasonable 
basis sets(Purvis and Bartlett, 1982; Stanton, 1997). Such 
methods are called ab initio^ but are not yet widespread for 
solids, where quantum Monte Carlo (QMC) is more often 
used(Filippi and Umrigar, 1996; Umrigar and Nightingale, 
1999). DFT calculations for molecules are usually much 
less computationally demanding, but the errors are less 
systematic and less reliable(Parr and Yang, 1995). 

However, many materials of current technological in¬ 
terest are both chemically complex and strongly corre- 
lated(Burke, 2012). Numerous metal oxide materials are 
relevant to novel energy technologies, such as Ti 02 for 
light-harvesting(0 ’Regan and Graetzel, 1991) or LiO com¬ 
pounds for batteries (Hafner et a7, 2011; Thackeray et a/., 
2012). For many cases, DFT calculations find ground- 
state structures and parameters, but some form of strong 
correlation method, such as introducing a Hubbard U or 
applying dynamical mean field theory (DMFT), is needed 
to correctly align bands and predict gaps (Anisimov et a7, 
1997a; Georges et a/., 1996). There is thus great in¬ 
terest in developing techniques that use insights from 
both ends, such as DFT+U and dynamical mean field 
theory (Anisimov et al.^ 1997b; Himmetoglu et a/., 2014; 
Kotliar et al.^ 2006; Kotliar and Vollhardt, 2004; Kulik 
and Marzari, 2010). 

There are two different approaches to combining DFT 
with lattice Hamiltonians (Capelle and Campo Jr., 2013). 
In the first, more commonly used, the lattice Hamilto¬ 
nian is taken as given, and a density function(al) the¬ 
ory is constructed for that Hamiltonian(Gunnarsson and 
Schonhammer, 1986). We say function(al), not functional, 
as the density is now given by a list of occupation num¬ 
bers, rather than a continuous function in real space. 
The parenthetical reminds us that although everything 
is a function, it is analogous to the functionals of real- 
space DFT. We will refer to this method as SOFT, i.e., 
site-occupation function(al) theory (Schonhammer et al, 
1995), although in the literature it is also known as lattice 
density functional theory (Ijas and Harju, 2010). While 
analogs of the basic theorems of real-space DFT can be 
proven such as the Hohenberg-Kohn (HK) theorems and 
the Levy constrained search formulation for SOFT, it is 
by no means clear (Harriman, 1986) how such schemes 
might converge to the real-space functionals as more and 
more orbitals (and hence parameters) are added. Alterna¬ 
tively, one may modify efficient solvers of lattice models 
so that they can be applied to real-space Hamiltonians 


(as least in 1-D), and use them to explore the nature of 
the exact functionals and the failures of present approxi- 
mations(Stoudenmire et a/., 2012; Wagner et a/., 2012). 
While originally formulated for Hubbard-type lattices, 
SOFT has been extended and applied to many different 
models include quantum-spin chains (Alcaraz and Capelle, 
2007), the Anderson impurity model(Carrascal and Ferrer, 
2012; Tows and Pastor, 2011), the 1-D random Fermi- 
Hubbard model (Xianlong et al.^ 2006b), and quantum 
dots(Schenk et al.^ 2011). 

These two approaches are almost orthogonal in philos¬ 
ophy. In the first, one finds approximate function(al)s for 
lattice Hamiltonians, and can then perform Kohn-Sham 
(KS) DFT calculations on much larger (and more inho¬ 
mogeneous) lattice problems (Campo and Capelle, 2005), 
but with all the usual caveats of DFT treatments (am 
I looking at interesting physics or a failure of an uncon¬ 
trolled approximation?). For smaller systems, one can 
often also compare approximate DFT calculations with 
exact results, results which would be prohibitively ex¬ 
pensive to calculate on real-space Hamiltonians. The 
dream of lattice models in DFT is that lessons we learn 
on the lattice can be applied to real-space calculations 
and functional developments. To this end, work has been 
done on understanding self-interaction corrections (Vieira 
and Capelle, 2010), and on wedding TDDFT and DMFT 
methods for application to more complex lattices (e.g. 
3-D Hubbard) (Karlsson et al.^ 2011a). And while it is 
beyond the scope of this current review, much work has 
been done on developing and applying density-matrix 
functional theory for the lattice as well(L6pez-Sandoval 
and Pastor, 2000, 2002, 2003, 2004; Saubanere and Pastor, 
2011, 2014). While such results can be very interesting, it 
is often unclear how failures of approximate lattice DFT 
calculations are related to failures of the standard DFT 
approximations in the real world. 

There is much interest in extracting excited-state infor¬ 
mation from DFT, and time-dependent (TD) DFT(Runge 
and Gross, 1984) has become a very popular first- 
principles approach (Burke et al.^ 2005; Marques et al.^ 
2012b; Ullrich, 2012). Because exact solutions and useful 
exact conditions are more difficult for TD problems, there 
has been considerable research using lattices. TD-SOFT 
can be proven for the lattice in much the same way SOFT 
is proven from ground-state DFT. This generalization is 
worked out carefully in Refs. (Farzanehpour and Tokatly, 
2012; Tokatly, 2011). Applications of TD-SOFT typically 
involve Hubbard chains both with and without various 
types of external potentials (Aryasetiawan and Gunnars- 
son, 2002; Karlsson et a/., 2011b; Mancini et al.^ 2014; 
Turkowski and Rahman, 2014). However, TD-SOFT has 
also been applied to the dimer to understand the effects of 
the adiabatic approximation in TD-DFT(Fuks et a/., 2013; 
Fuks and Maitra, 2014a,b), strong correlation (Turkowski 
and Rahman, 2014), and TD-LDA results for stretched 
H 2 in real-space (Aryasetiawan et a/., 2002). Unfortu- 


3 


nately, we will already fill this article simply discussing 
the ground-state SOFT problem, and save the TD case 
for future work. 
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FIG. 1 Many-body view of two distinct regimes of the asym¬ 
metric Hubbard dimer. On the left, the charging energy is 
much greater than the difference in on-site potentials. On the 
right, the situation is reversed. 

To get the basic idea, consider Fig. 1. It shows the 
asymmetric Hubbard dimer in two different regimes. On 
the left, the Hubbard U energy is considerably larger than 
the difference in on-site potentials and the hopping energy 
t. This is the case most often analyzed, where strong cor¬ 
relations drive the system into the Mott-Hubbard regime 
if U is also considerably larger than t. The on-site occu¬ 
pations are in this case close to 1. On the right panel, 
U is in contrast smaller than the on-site potential differ¬ 
ence and here the dimer stays in the charge-transfer 
regime, where both electrons mostly sit in the same deeper 
well. This is the many-body view of the physics of an 
asymmetric Hubbard dimer. 
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FIG. 2 Occupations n and potentials v of an asymmetric 
half-filled Hubbard dimer as a function of U. The on-site 
potential difference Av is shown in black and the KS on-site 
potential difference Ai’s is in red. The second and third panels 
correspond to the situations of Fig. 1. 

Now we turn to the KS-DFT viewpoint. Here, we 


replace the interacting Hubbard dimer {U ^ 0) with a 
non-interacting {U = 0) tight-binding dimer, called the 
KS system, that reproduces the Hubbard occupations. 
In Fig. 2, we take the asymmetric dimer with the same 
on-site potential difference, but we vary U. We plot the 
occupations, showing how, as U increases, their difference 
decreases. But we also plot the on-site potentials of the 
Kohn-Sham model, Avq^ that are chosen to reproduce the 
occupations of the interacting system with a given value 
of U. As U increases, the KS on-site potential difference 
reduces and the offset from 0 increases. The middle panel 
corresponds to the charge-transfer conditions of Fig. 1, 
while the last panel corresponds to the Mott-Hubbard con¬ 
ditions of Fig. 1. The basic theorems of DFT show that 
if we know the energy as a function(al) of the density, we 
can determine the occupations by solving effective tight- 
binding equations, the KS equations, and then find the 
exact ground-state energy. This is not mean-field theory. 
It is instead a horribly contorted logical construction, that 
is wonderfully practical for computations of ground-state 
quantities. Inside this article, we give explicit formulas 
for the energy functional of the Hubbard dimer. 

We perform a careful study of the Hubbard dimer, to 
show the differences between SOFT and real-space DFT. 
We show how it is necessary to introduce inhomogeneity 
into the site occupations in order to find the exact density 
function(al) explicitly. In Section 2.A we explain the logic 
of the KS DFT approach in excruciating detail in order 
to both illustrate the concepts to those unfamiliar with 
the method and to give explicit formulas for anyone doing 
SOFT calculations. We elucidate the differences between 
the KS and the many-body Green’s functions in Section 
4.C. Next, in Sections 4 and 5 we discuss in detail both 
concepts and tools for strong correlation, and explain 
how the gap problem appears in DFT. We construct the 
adiabatic connection formula for the exact function(al) in 
Section 5.B, showing how it is quantitatively similar to 
those of real-space DFT. We use the theory to construct 
a simple parameterization for the exact function(al) for 
this problem in Section 6, where we also demonstrate 
the accuracy of our formula by finding ground-state ener¬ 
gies and densities by solving the KS equations with our 
parametrization. In Section 7.A, we study the broken- 
symmetry solutions of Hartree-Fock theory, showing that 
these correctly yield both the strongly-correlated limit 
and the approach to this limit for strong correlation. In 
Section 7.B we present BALDA (Bethe-ansatz local den¬ 
sity approximation), a popular approximation for lattice 
DFT, and in Section 7.C we compare the accuracy of 
BALDA and Hartree-Fock to each other. We discuss 
fractional particle number and the derivative discontinu¬ 
ity in Section 8. Finally, we end with a discussion of 
our results in Section 9. In Table I we list our notation 
for the Hubbard dimer, as well as many standard DFT 
definitions. 

Our purpose here is several-fold. Perhaps most im- 
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Definition 

Description 

Generic DFT 

T[n] 

Many-body wfn of density n 

4>[n] 

Kohn-Sham wfn of density n 

F = T+Vee 

Hohenberg-Kohn Functional 

Exc = F - Ts - Uh 

Exchange-correlation energy 

Ex = - Uh 

Exchange energy 

Ex = -Un/2 

Exchange energy for 2 electrons 

Ec = Tc Uc 

Total correlation energy 

Tc=T-Ts 

Kinetic correlation energy 

Uc = Ke -Un-Ex 

Potential correlation energy 

Uxc(X) = t/^c/A 

Adiabatic connection integrand 

Tc = Ec-dE^/dX\x=i 

Method to extract Tc from Ec 

Uc = dE^/dX\x=i 

Method to extract Uc from Ec 

hs — — 

Kohn-Sham hamiltonian 

i;s = -c + 'Ch + 'Cxc 

Kohn-Sham one-body potential 

^trad _ ^HF 

Quantum chemical corr. energy 

SOFT Hubbard 

ni, 712 

Occupations at sites 1, 2 

iV = ni + 772 

Total number of electrons 

(M 

1 

II 

<1 

Occupation difference 

AtTI = 7771 — r772 

Magnetization difference 

V\, V2 

On-site potentials 

V — V2)/2 — 0 

On-site potential average 

At; = t;2 — t;i 

On-site potential difference 

At;xc = 'Cxc ,2 - 'Cxc,i 

XC potential difference 

Un = U{N‘^ + An^)/A 

Hartree energy 

Fhx = U{N‘^ + A772)/8 

Hartree-Exchange energy 

Ts =-V(2-|iV-2|)2- 

-A77^ Single particle hopping energy 

Dimensionless Variables 

t = E/2t 

Energy in units of hopping 

u = U/2t 

Hubbard U in units of hopping 

V = Av/2 1 

Pot. diff. in units of hopping 

p = An/2 

Reduced density difference 

p=\- p 

Asymmetry parameter 


TABLE I Standard DFT definitions and our Hubbard dimer 
notation. 

portantly, this article is intended to explain the logic of 
modern DFT to our friends who are more familiar with 
strongly correlated systems. We take the simplest model 
of strong correlation, and illustrate many of the basic 
techniques of modern DFT. There are many more tricks 
and constructions, but we save those for future work. The 
article should be equally useful to researchers in other 
fields who are unfamiliar with the logic of DFT, such as 
traditional quantum chemists or atomic and molecular 
physicists. 

Secondly, the article forms an essential reference for 
those researchers interested in SOFT, possibly in very 
different contexts and applied to very different models. It 
shows precisely how concepts from first-principles calcula¬ 
tions are realized in lattice models. Third, we give many 
exact results for this simple model, expanding in many dif¬ 
ferent limits, showing that even in this simple case, there 
are orders-of-limits issues. Fourth, we use DFT techniques 
to find a simple but extremely accurate parametrization 
of the exact function(al) for this model. Even though the 


model can be solved analytically, the function(al) cannot 
be expressed explicitly. Thus our parametrization pro¬ 
vides an ultra-convenient and ultra-accurate expression 
for the exact function(al) for this model, that can be used 
in the ever increasing applications of SOFT. Finally, we 
examine several standard approximations to SOFT, in¬ 
cluding both restricted and unrestricted mean field theory, 
and the BALDA, and we find surprising results. 

2. BACKGROUND 

In this section we briefly introduce real-space DFT, 
and the logical underpinnings for everything that follows. 
Then we discuss the mean-field approach to the Hubbard 
model as well as a few well-known results and limits for the 
Hubbard dimer. Throughout this section we use atomic 
units for all real-space expressions so all energies are in 
Hartree and all distances are in Bohr. 

A. Density functional theory 

We restrict ourselves to non-relativistic systems within 
the Born-Oppenheimer approximation with collinear mag¬ 
netic fields(Engel and Dreizler, 2011). Density functional 
theory is concerned with efficient methods for finding the 
ground-state energy and density of N electrons whose 
Hamiltonian contains three contributions: 

^ = T + We + (1) 

The first of these is the kinetic energy operator, the second 
is the electron-electron repulsion, while the last is the one- 
body potential, 

N 

C = ^t>(ri). (2) 

Only N and u(r) change from one system to another, be 
they atoms, molecules or solids. In 1964, Hohenberg and 
Kohn proved that for a given electron-electron interaction, 
there was at most one u(r) that could give rise to the 
ground-state one-particle density no(r) of the system, 
thereby showing that all ground-state properties of that 
system were uniquely determined by no(r) (Hohenberg 
and Kohn, 1964). The ground-state energy Eq could then 
be found by splitting the variational principle into two 
steps via the Levy-Lieb constrained search approach (Levy, 
1979; Lieb, 1983). First, the universal functional F is 
determined, 

F[n] = min(^| f + We |^) = T[n] + WeW (3) 

where the minimization is over all normalized, antisym¬ 
metric T with one-particle density n(r). This establishes a 
one-to-one connection between wavefunctions and ground- 
state densities, and enables us to define the minimizing 
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wavefunction functional ^[no]. Then the ground-state 
energy is determined by a second minimization step of 
the energy functional E[n]^ 


+ J d^rn{r)v{r) 

This shows that Eq can be found from a search over 
one-particle densities n(r) instead of many-body wave- 
functions T, provided that the functional E[n] is known. 
The Euler equation corresponding to the above minimiza¬ 
tion for fixed N is simply 



Eq = min {E[n\} = min < E[n\ 


6E[n] 

6n{r) 


no(r) 


-v{r). 


( 5 ) 


Armed with the exact T[n], the solution of this equation 
yields the exact ground-state density which, when inserted 
back into F[n], yields the exact ground-state energy. 

To increase accuracy and construct E[n], modern DFT 
calculations use the Kohn-Sham (KS) scheme that imag¬ 
ines a fictitious set of non-interacting electrons with the 
same ground-state density as the real Hamiltonian(Kohn 
and Sham, 1965). These electrons satisfy the KS equa¬ 
tions: 


+ Vs(r)| = €i (6) 

where Vs{r) is defined as the unique potential that gen¬ 
erates single-electron orbitals 0i(r) that reproduce the 
ground-state density of the real system. 


^o(r) = y] |</>i(r)p. (7) 

occ 

To relate these to the interacting system, we write 


E[n\ — Ts[n] -^Uu[ti] + ExcM- (^) 

Tg is the non-interacting (or KS) kinetic energy, given by 

1 r ^ 

Ts{n] = 2 / ( 9 ) 


where we have assumed the KS wavefunction (as is almost 
always the case) is a single Slater determinant ^ of single¬ 
electron orbitals. The second expression follows from Eq. 
(3) applied to the KS system, it emphasizes that Tg is a 
functional of n(r), and the minimizer defines ^[no], the 
KS wavefunction as a density functional. Then 1 /h[^] is 
the classical electrostatic self-repulsion of n(r). 


Unln] = j d- 


,3 , n(r) n{r') 
Ir-r'| ^ 


( 10 ) 


and Exc is called the exchange-correlation energy, and is 
defined hy Eq. (8). 


Lastly, we differentiate Eq. (8) with respect to the 
density. Applying Eq. (5) to the KS system tells us 


Vs{r) = - 


<5rs[n] 
6n{r) ’ 


(11) 


yielding 

Vs{r) = v{r) + i;H(r) + i;xc(r) (12) 

where VH(r) is the classical electrostatic potential and 


t’xc(r) 


SExc 

6n{r) 


(13) 


is the exchange-correlation potential. This is the single 
most important result in DFT, as it closes the set of 
KS equations. Given any expression for Exc in terms of 
no(r), either approximate or exact, the KS equations can 
be solved self-consistently to find no(r) for a given i;(r). 

However, we also note that, just as in all such schemes, 
the energy of the KS electrons does not match that of the 
real system. This ‘KS energy’ is 


Es[n]=y2^i = Ts + Vs, (14) 

i 


but the actual energy is 

Eo = E[no] + V[no] = Ts[no] + Uu[no] + E^c[no] + V[no] 

(15) 

where no(r) and Ts[no] have been found by solving the 
KS equations, and inserted into this expression. Thus, in 
terms of the KS orbital energies, there are double-counting 
corrections, which can be deduced from Eqs. (14) and 
(15): 

Eo = Es- Uii[no] EE^c[no] - / d^r no{r) v^c[no]{r). 

(16) 

We emphasize that, with the exact E^c[i^o]^ solution of 
the KS equations yields the exact ground-state density 
and energy, and this has been done explicitly in model 
cases(Wagner et al.^ 2013), but is computationally ex¬ 
orbitant. The practical use of the KS scheme is that 
simple, physically motivated approximations to Exc[(^o] 
often yield usefully accurate results for Eq, bypassing 
direct solution of the many-electron problem. 

For the remainder of this article, we drop the subscript 
0 for notational convenience, and energies will be assumed 
to be ground-state energies, unless otherwise noted. For 
many purposes, it is convenient to split E^c into a sum 
of exchange and correlation contributions. The definition 
of the KS exchange energy is simply 

Ex [n] = ($ [n] I fee I $ [n] )-U^[n], (17) 

The remainder is the correlation energy functional 

Ea[n] = F[n]- ($[n]| f + Ke |$[n]), (18) 
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which can be decomposed into kinetic and potential 
Uc contributions (see Eqs. (76) and (77) in Sec. 5). Addi¬ 
tionally, all practical calculations generalize the preceding 
formulas for arbitrary spin using spin-DFT (von Barth 
and Hedin, 1972). 

For just one particle (N = 1), there is no electron- 
electron repulsion, i.e., Fee = 0- This means 

Fx = -[/h, Ec = 0, {N = 1), (19) 

i.e., the self-exchange energy exactly cancels the Hartree 
self-repulsion. Since there is no interaction, F^[n] = 
T[n] = Ts[n], and for one electron we know the explicit 
functional: 

Ts = = J (fr |Vn| V(8n), (20) 

which is called the von Weisacker functional (Weizsacker, 
1935). For two electrons in a singlet {N = 2), 

E^ = -U^I2, Ts = tW, (iV = 2), (21) 

but the correlation components are non-zero and non¬ 
trivial. 

Many popular forms of approximation exist for E^c[n]^ 
the most common being the local density approxima¬ 
tion (LDA)(von Barth and Hedin, 1972; Kohn and Sham, 
1965; Perdew and Wang, 1992), the generalized gradi¬ 
ent approximation (GGA)(Becke, 1988; J. P. Perdew and 
Fiolhais, 1992; Lee et al.^ 1988; Perdew, 1986a; Perdew 
et al.^ 1996a), and hybrids of GGA with exact exchange 
from a Hartree-Fock calculation (Adamo and Barone, 1999; 
Becke, 1993a; Heyd et a/., 2003; Perdew et a/., 1996b). 
The computational ease of DFT calculations relative to 
more accurate wavefunction methods usually allows much 
larger systems to be calculated, leading to DFT’s immense 
popularity today (Pribram-Jones et a/., 2015). However, 
all these approximations fail in the paradigm case of 
stretched 772, the simplest example of a strongly corre¬ 
lated system(Baerends, 2001; Hellgren et a/., 2014). 

B. The Hubbard model 

The Hubbard Hamiltonian is possibly the most stud¬ 
ied, and simplest, model of a strongly correlated elec¬ 
tron system. It was initially introduced to describe the 
electronic properties of narrow-band metals, whose con¬ 
duction bands are formed by d and / orbitals, so that 
electronic correlations become important (Fradkin, 2013; 
Hubbard, 1963). The model was used to describe fer¬ 
romagnetic, antiferromagnetic and spin-spiral instabili¬ 
ties and phases, as well as the metal-insulator transition 
in metals and oxides, including high-Tc superconduc- 
tors(Dagotto, 1994; Lee et a/., 2006). The Hubbard model 
is both a qualitative version of a physical system depend¬ 
ing on what terms are built in( Anderson, 1987; Schulz, 


1990) and also a testing-ground for new techniques since 
the simpler forms of the Hubbard model are understood 
very well(Bickers and Scalapino, 1989; Bickers et a/., 1989; 
Hirsch, 1989, 1993). 

The model assumes that each atom in the lattice has 
a single orbital. The Hamiltonian is typically written as 
(Essler et a/., 2005; Ha, 1996; Mattis, 1993; Takahashi, 
2005) 

H ^ ^ Via ^ ^ a ^ T T ^ ^ Ui fli'^ 

i,cr i j (J i 

( 22 ) 

where at its simplest the on-site energies are all equal 
Via- = 0 as well as the Goulomb integrals Ui = U. Further, 
the hopping integrals tij typically couple only nearest 
neighbor atoms and are equal to a single value t. 

We note that here the interaction is of ultra-short range, 
so that two electrons only interact if they are on the same 
lattice site. Further, they must have opposite spins to 
obey the Pauli principle. Simple examples of building 
in more complicated physics include using next-nearest- 
neighbor hoppings or nearest neighbors Goulomb integrals 
for high-Tc cuprate calculations and magnetic proper- 
ties(Daul and Noack, 1998; Duffy and Moreo, 1995; Lin 
and Hirsch, 1987), and varying on-site potentials used 
to model confining potentials (Robert et a/., 2008). Also, 
adding more orbitals per site delivers multi-band Hub¬ 
bard models, where Coulomb correlations may be added 
to some or all of the orbitals. The Hubbard model has 
an analytical solution in one dimension, via Bethe ansatz 
techniques (Lieb and Wu, 2003, 1968). 

If the Hubbard U is small enough, a paramagnetic mean- 
field (MF) solution provides a reasonable description of 
the model in dimensions equal or higher than two. As an 
example, the Hubbard model in a honeycomb lattice can 
describe correctly a number of features of gated graphene 
samples (Herbut, 2006). However, for large U or in one 
dimension, more sophisticated approaches are demanded, 
which go beyond the scope of this article(Fradkin, 2013; 
Lieb and Wu, 1968). 

We describe briefly the well-known broken-symmetry 
MF solution, where the populations of up- and down-spin 
electrons can differ. The standard starting point for the 
MF solution neglects completely quantum fluctuations: 

(ni-|- - n,t) (fiii - riii) = 0, (MF) (23) 

where ni<j = (fiia), so that 

kf ^ = F ^ '^a) ■ (24) 

i 

The MF hamiltonian is then just an effective single¬ 
particle problem 

= (25) 

i<j 

nia - t F(Tcia + h.c.), (26) 

j 
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where = Vi^ -\- U rii^. This can be easily 

diagonalized if one assumes space-homogeneity of the 
occupations rii^a- = For large [/, the broken symmetry 
solution (often ferromagnetic) has lower energy than the 
paramagnetic solution. 


C. The two-site Hubbard model 


We now specialize to a simple Hubbard dimer model 
with open boundaries, but we allow different on-site spin- 
independent energies by introducing a third term that 
produces asymmetric occupations, 

H = {c\^C 2 cr + h.c) + u'^hif hii + y] Vihi (27) 

cr i i 


where we have made the choices ti 2 = 1 21 = t and U 1 +U 2 = 
0. Our notation for this Hamiltonian can be found in 
Table 1. 



FIG. 3 Ground-state energy of Hubbard dimer as a function 
of Av for several values of U and 2t = 1. 


It is straightforward to find an analytic solution of the 
model for any integer occupation N. However, we special¬ 
ize to the particle sub-space N = 2, Sz = 0 in what follows 
unless otherwise stated. We expand the Hamiltonian in 
the basis set [|1 t 1 I}, |11 2 \1 i 2 t}, |2 t 2 ;}]: 


( 2v\ -^-U —t t 0 \ 

-t 0 0 -t 

t 0 0 t 

y 0 —tt 2 v2 + U j 


(28) 


The eigenstates are three singlets and a triplet state. 
The ground-state energy corresponds to the lowest-energy 
singlet, and can be found analytically. The expressions are 
given in appendix A. The wavefunction, density difference, 
and individual energy components are also given there. 
We plot in Fig. 3 the ground-state energy as a function 
of Av for several values of t/, while in Fig. 4, we plot the 
occupations. 



FIG. 4 Ground-state occupation of Hubbard dimer as a func¬ 
tion of Av for several values of U and 2t = 1. 

When U = 0, we have the simple tight-binding result, 
for which the ground-state energy is 

E = + Av^ iU = 0), (29) 

An = 2 An/v/(2t)2 + {U = 0). (30) 

where An is defined in Table I. If there is only one electron, 
these become smaller by a factor of 2. The curves for U = 
0.2 are indistinguishable (by eye) from the tight-binding 
result. We may simplify the expressions by introducing 
an effective hopping parameter, 

i = t^yl + {Av/{2t))^ (31) 

which accounts for the asymmetric potential. Then 

E = -2i {U = 0), (32) 

An = Av/i {U = 0), (33) 

i.e., the same equations as when Av = 0. 

In the other extreme, as U grows, we approach the 
strongly correlated limit. For a given Au, as U increases. 
An decreases as in Figs. 2 and 4, and the magnitude of 
the energy shrinks. Typically, the E[Av) curve morphs 
from the tight-binding result towards two straight lines 
for U large: 

E-{U- Av) e{Av -U) U'>2t (34) 

Anc^2e{Av -U) [/>2t (35) 

We also have a simple well-known result for the symmetric 
limit, Au=0, where 

E = -^/{2i^V(UjW+Ul2 {An = Av = 0) (36) 

This vanishes rapidly with 1/U for large U. Its behavior is 
different from the case with finite Av. Results for various 
limits and energy components are given in Appendix A. 


D. Quantum chemistry 

Traditional quantum chemical methods (often referred 
to as ab initio by their adherents) usually begin with 







































the solution of the Hartree-Fock equations(Szabo and 
Ostlund, 1996). For our Hubbard dimer, these are noth¬ 
ing but the mean-field equations of Sec 2.B. Expressing 
the paramagnetic HF Hamiltonian of Eq. (26) for two 
sites yields a simple tight-binding Hamiltonian and eigen¬ 
value equation describing a single-particle in an effective 
potential: 

v'f{ni) = Vi +Uni/2. (37) 

with an eigenvalue: 

=(u- ^(Av«ff)2 + (2i)2^ /2. (38) 

Writing 0®^ = (ci,C 2 )^, then 

An = 2 (cf — c|) = 2 ^ ^ 

where x = Av^^/2t, and ^ + 1 — Eq. (39) is 

quartic in An and can be solved algebraically to find An 
as a function of Av explicitly (appendix E). Just as in 
KS, the HF energy is not simply twice the orbital energy, 
there is a double-counting correction: 

^MF ^ - Un (40) 



These energies are plotted in Fig 5. We see that for 



FIG. 5 Ground-state energy of the Hartree-Fock Hubbard 
dimer (thick dashed line) and exact ground-state of the Hub¬ 
bard dimer (thin solid line) as a function of Av for several 
values of U and 21 = 1. 

small U, HF is very accurate, but much less so for 21 ^ 
U <C Av. In fact, the HF energy becomes positive in 
this region, unlike the exact energy, which we prove is 
never positive in appendix C. The molecular orbitals 
often used in chemical descriptions have traditionally 
been those of HF calculations, despite the fact that HF 
energies are usually far too inaccurate for most chemical 


energetics (Bickelhaupt and Baerends, 2000). (They have 
now largely been supplanted by KS orbitals.) In quantum 
chemical language, the paramagnetic mean-field solution 
is called restricted HF (RHF) because the spin symmetry 
is restricted to that of the exact solution, i.e., Sz = 0. For 
large enough t/, the broken-symmetry, or unrestricted, 
solution is lower, and is labeled UHF, which we discuss 
in Sec. 7.A. 



Av 

FIG. 6 Gorrelation energy of Hubbard dimer as a func¬ 

tion of Av for several values of U and 2t = I. 

Accurate ground-state energies, especially as a function 
of nuclear positions, are central quantities in chemical elec¬ 
tronic structure calculations(Szabo and Ostlund, 1996). 
Most such systems are weakly correlated unless the bonds 
are stretched. The correlation energy of traditional quan¬ 
tum chemistry is defined as just the error made by the 
(restricted) HF solution: 

^trad ^ ^ _ ^HF_ 

This is plotted in Fig 6. This is always negative, by the 
variational principle. Many techniques have been highly 
developed over the decades to go beyond HF. These are 
called model chemistries, and for many small molecules, 
errors in energy differences of less than I kcal/mol (0.05 
eV) are now routine (Bartlett and Musial, 2007; Ochterski 
et a/., 1995). 

Usually is a small fraction of E for weakly corre¬ 

lated systems. For example, for the He atom, E = —77.5 
eV, but E^^^ = —1.143 eV. This is the error made by 
a HF calculation. In Fig. 6 we plot E^^^ just as we 
plotted E in Fig. 5. We see that for strong correlation 
E^^^ becomes large (^ —Uf2 for Av ^ U), much larger 
than E. However, E is much smaller, and so any strongly 
correlated method should reproduce E accurately. In fact, 
one can already see difficulties for weakly correlated ap¬ 
proximations in this limit. For weak correlation, a small 
percent error in E^^^ yields a very small error in E, but 
produces an enormous error in E in the strong correlation 
limit. For an infinitely stretched molecular bond, t —> 0 
while U remains finite, so only one electron is on each site. 
Thus E ^ 0, so we can think of E as the ground-state 
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electronic energy relative to the dissociated limit, i.e. the 
binding energy. 

Because HF is accurate when correlation is weak, and 
because quantum chemistry focuses on energy differences, 
the error is often measured in terms of the accuracy of the 
exchange-correlation together (if both are approximated 
as in most DFT calculations). For 2 electrons having 
5';^ = 0, the exact exchange is trivial, and so we will focus 
on approximations to the correlation energy. 

Notice the slight difference in definition of correlation 
energy between DFT (Eq. 18) and quantum chemistry 
(Eq. (41)) (Gross et al.^ 1996; Sahni and Levy, 1986; 
Umrigar and Gonze, 1994). In DET, all quantities are 
defined on a given density, usually the exact density of the 
problem, whereas in quantum chemistry, the HE energy 
is evaluated on the density that minimizes the HF energy. 
For weakly correlated systems, this difference is extremely 
small(Gorling and Ernzerhof, 1995), but is not so small 
for large U. And, one can prove, (Gross 

et a/., 1996), (see appendix C). 

We close by emphasizing the crucial difference in phi¬ 
losophy between DFT and traditional approaches. In 
many-body theory, mean-field theory is an approximation 
to the many-body problem, yielding an approximate wave- 
function and energy which are expected to be reasonably 
accurate for small U. In DFT, this treatment arises from 
approximating F for small I/, and so should yield an ac¬ 
curate KS wavefunction and expectation values for small 
U. Thus, only one-body properties that depend only on 
position are expected to be accurate, and their accuracy 
can be improved by further improving the approximation 
to F. For large t/, such an approximation fails, but there 
is still an exact F that yields an exact answer. 


3. SITE-OCCUPATION FUNCTION THEORY (SOFT) 

In this section, we introduce the site-occupation func¬ 
tion theory for the Hubbard dimer (Carrascal and Ferrer, 
2012; Gunnarsson and Schonhammer, 1986; Schonhammer 
and Gunnarsson, 1988; Schonhammer et a/., 1995). If we 
want a physical system where this arises, think of stretched 
H 2 (McLean et a/., 1960). We imagine a minimal basis set 
of one function per atom for the real Hamiltonian. We 
choose these basis functions to be 1 5 orbitals centered on 
each nucleus, but symmetrically orthonormalized. Then 
each operator in real-space contributes to the parameters 
in the Hubbard Hamiltonian as seen in Appendix F. 

It is reasonably straightforward to establish the validity 
of SOFT for our dimer. So long as each occupation can 
come from only one value of Ai;, for a fixed 1/, there is a 
one-to-one correspondence between An and Ai;, and all 
the usual logic of DFT follows. But note that T and V in 
SOFT do not correspond to the real-space kinetic energy 
and potential energy. For example, the hopping energy is 
negative, whereas the real-space kinetic energy is positive. 


This means that all theorems of DFT to be used must 
be reproven for the lattice model. More importantly, the 
SOFT does not become real-space DFT in some limit of 
complete basis sets (in any obvious way). We will however 
apply the same logic as real-space DFT, with the hopping 
energy in SOFT playing the role of the kinetic energy in 
DFT, and the on-site energy in SOFT playing the role of 
the one-body potential. The interaction term obviously 
plays the role of 


A. Non-interacting warm-up exercise 


To show how SOFT works, begin with the U = 0 case, 
i.e., tight-binding of two non-interacting electrons. The 
ground-state is always a spin singlet. From the non¬ 
interacting solution, we can solve for Av in terms of An 


2 1 An 
Vd — An^ ’ 


(42) 


and substitute back into the kinetic energy expectation 
value to find 


T(ni,n2) = -2ty/nin2. (43) 

This is the universal density function(al) for this non¬ 
interacting problem (see Eq. (3)), and can be used to 
solve every non-interacting dimer. 

To solve this N = 2 problem in the DFT way, we note 
that T is playing the role of F( 77 . 1 , 722 ). So the exact 
function(al) here is 

F{ni) = -2ty/nin2, {U = 0) (44) 

from which we can calculate all the quantities of interest 
using a DFT treatment. Note that everything is simply a 
function(al) of 77.1 since 77.2 = {N — 77 . 1 ), or alternatively a 
function(al) of An. 

We then construct the total energy function(al): 

E{ni) = F{ni) — Av An/2^ {U = 0) (45) 

and minimize with respect to 77.1 for a given Av to find 
the ground-state energy and density: 

E = - v/(2i)2 + Av2, (46) 

An = 2 Av/ v/( 2 i )2 + A?;2. ( 47 ) 

Both of these agree with the traditional approach and 
recover Eqs. (29) and (30). The N = 1 result is half as 
great as Eqs. (46) and (47). 

We can deduce several important lessons from this 
example. Eirst, we need to vary the one-body potential 
(in this case, the on-site energy difference) to make the 
density change through all possible values, in order to 
find the function(al), since it requires knowing the one- 
to-one correspondence for all possible densities. Second, 








10 


if we really change the atoms in our 2-electron stretched 
molecule, of course the minimal basis functions would 
change, and both t and Av would differ. But here we keep t 
fixed, and vary Av simply to explore the function(al), even 
if we are only interested in solving the symmetric problem. 
(Real-space DFT does not suffer from this problem, as 
the kinetic and repulsion operators are universal.) Third, 
we are reminded that the hopping and on-site operators 
in no sense represent the actual kinetic and one-body 
potential terms - they are a mixture of each. Finally, 
although we ‘cheated’ and extracted the kinetic energy 
function(al) from knowing the solutions, if someone had 
given us the formula, it would allow us to solve every 
possible non-interacting Hubbard dimer by minimizing 
over densities. And an approximation to that formula 
would yield approximate solutions to all those problems. 


B. The interacting functional 

For the interacting case, we cannot analytically write 
down the exact function(al) F{ni) at = 2 in closed 
form. Although we have analytic formulas for both E 
and An as functions of Av^ the latter cannot be explicitly 
inverted to yield an analytic formula for F{An). However, 
we can plot the function(al), by simply plotting F = E — V 
as a function of ni, and see how it evolves from the U = 0 
case to stronger interaction. The spin state is always a 
singlet. We plot in Fig. 7 the F-function(al) as a function 



0 0.5 1 1.5 2 


FIG. 7 F-function(al) of Hubbard dimer as a function of ni 
for several values of U and 21 = 1. 


The oldest form of DFT (Thomas-Fermi theory(Fermi, 
1928; Thomas, 1927)) approximates both T{ni) and 
Hee(^i) and so leads to a crude treatment of the en¬ 
ergetics of the system. A variation on this was used in 
Ref. (Campo and Capelle, 2005) to enable extremely 
large calculations. 


C. Kohn-Sham method 

The modern world uses the KS scheme, and not pure 
DFT(Burke, 2012). The scheme in principle allows one 
to find the exact ground-state energy and density of an 
interacting problem by solving a non-interacting one. This 
scheme is what produces such high accuracy while using 
simple approximations in DFT calculations today. Next, 
we see how the usual definitions of KS-DFT should be 
made for our dimer. 

The heart of the KS method is the fictitious system 
of non-interacting electrons whose density matches with 
the ground-state density of the interacting system. For 
our two-electron system, the KS system is that of non¬ 
interacting electrons {U = 0) with an on-site potential 
difference Avq^ defined to reproduce the exact An of the 
real system. This is just the tight-binding problem with 
an effective on-site potential difference, and is illustrated 
in Fig. 2. 

As stated in Section 2.A, in KS-DFT one convention¬ 
ally extracts the Hartree contribution from the electron- 
electron repulsion. There are deep reasons for doing 
so, which center on the remnant, the XC energy, be¬ 
ing amenable to local and semilocal-type approxima- 
tions(Burke et a/., 1998; Pribram-Jones et a/., 2015). To 
see how the Hartree energy should be defined here, rewrite 
the electron-electron repulsion as: 

G = ^ ~ ^4)- (^0) 

i 

This form mimics the treatment in DFT. The first term 
depends only on the total (i.e. spin-summed) density, akin 
to Hartree in real-space DFT. The remaining terms cancel 
the self-interaction that arises from using the total density 
for the electron-electron interaction. For the N = 2 dimer, 
this decomposition results in 

f7H(An) = ^ (nf + nl) , (51) 


of ni for several values oiU. As U increases we can see 
F appears to tend to U\1 — nfi. 

For any real problem the Euler equation for a given Av 
is 


dF{ni) Av 
dni 2 


(48) 


and 

E^{An) = {nf + nj) , (52) 

which satisfies = —Unl2> for N = 2 as defined in real- 
space DFT for a spin singlet, Eq. (23). Together, the 
Hartree-Exchange is 


and the unique ni{Av) is found that satisfies this. Then 
E{Av) = F{ni, Av) — AvAn{Av)/2. (49) 


E^^{An) = j{nl + nl) = j . 


(53) 
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In Appendix B we see that the leading order in the U 
expansion of the F—function(al) yields the same result. 
A typical mean field treatment of also results in Eq. 
(53). In DFT there is always self-exchange, even for 
one or two particles. In many-body theory, exchange 
means only exchange between different electrons. Despite 
this semantic difference, both approaches yield the same 
leading-order-in-expression for the dimer, which we call 
Fhx here (but is often called just Hartree in many-body 
theory). 

For the dimer, from Eq. (43), the KS kinetic energy is 
just 

Ts{ni) = -2ty/nin2, (54) 

so that F^^(ni) = Ts(ni) -h Fhx(^i) as in Section 2.D. 
We can then define the correlation energy function from 
Eq. (18), so that 

Fc(ni) = F(ni) - Ts(ni) - (55) 

In Fig. 8, we plot the correlation energy as a function of 



FIG. 8 Plot of exact Ec (blue line) and Fc,par (red dashed 
line) for different U and 21 = 1. 

rii. For small I/, 

U <^2t (56) 

which is much smaller than the Hartree-exchange contri¬ 
bution, and is a relatively small contribution to E. But 
as U increases, 

Ec - -U{1 - (ni - 1)2)/2, f7 » 21 (57) 

with a cusp at half-filling. Combined with Fhx, this 
creates E for large U as in Fig. 7. 

Inserting this result into Eq. (48), we find that the KS 
electrons have a non-interacting Hamiltonian: 

hs\4>) = <^s\4>), (58) 


where this KS Hamiltonian is 

hs{An) = -t (c|c 2 + ft..c.) + y] Va,i{An)hi. (59) 

i 

The KS potential difference is 

Aus(An) = Au — I/An/2 + Auc(An), (60) 

where 

Auc = —2dFc(ni)/dni. (61) 

This is the key formal result of the KS paper (Kohn and 
Sham, 1965) as applied to SOFT: For any given form of 
the (exchange-)correlation energy, differentiation yields 
the corresponding KS potential. If the exact expression 
for Fc(ni) is used, this potential is guaranteed (Wagner 
et a/, 2013) to yield the exact ground-state density when 
the KS equations are iterated to convergence via a simple 
algorithm. 

In Fig. 9, we plot several examples of the dependence 
of the potentials in the KS system as a function of ni, 
which range from weakly (JJ = 0.4) to strongly {U = 
10) correlated cases. In each curve, the black line is 
the actual on-site potential difference as a function of 
occupation of the first site. The blue line is the KS 
potential difference, which is the on-site potential needed 
for two non-interacting (1/ = 0) particles to produce the 
given ni. This is found by inverting the tight-binding 
equation for the density, Eq. (42). Their difference is the 
Hartree-exchange-correlation on-site potential, denoted 
by the red line. Finally, the green line is just Hartree- 
exchange, which ignores correlation effects. For U = 0.4, 
we see that the difference between blue and black is quite 
small, and almost linear. Indeed the Hartree-exchange 
contribution is always linear (see Eq. (60)). Here the 
red is indistinguishable by eye from the green, showing 
how small the correlation contribution to the potential is. 
This means the HF and exact densities will be virtually 
(but not quite) identical. When we increase U to I, we 
see a similar pattern, but now the red line is noticeably 
distinct from the green. For any given ni, the blue curve 
is smaller in magnitude than the black. This is because 
turning on U pushes the two occupation numbers closer, 
and so their KS on-site potential difference is smaller. 
Again, the red curve is larger in magnitude than the 
green, showing that HF does not suppress the density 
difference quite enough. In our final panel, U = 10, and 
the effects of strong correlation are clear. Now there is a 
huge difference between black and blue curves. Because 
U is so strong, the density difference is close to zero for 
most ni, making the blue curve almost ffat except at the 
edges. In the KS scheme, this is achieved by the red curve 
being almost ffat, except for a sudden change of sign near 
ni = I. These effects give rise to the Avq values shown 
in Fig. 2. This effect is completely missed in HF. 
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FIG. 9 Plots of Avs (blue) and its components, Av (black), 
—f/An/2 (green), and Avc — UAnj2 (red) plotted against ni 
for various U and 2t — 1. The arrows indicate the occupations 
used in Fig. 2. 

To emphasize the role of correlation, in Fig. 10, we plot 
the correlation potential alone, which is the difference 
between the red and green curves in Fig. 9. Values from 
the blue curves for Av = 2 were used to make Fig. 2. Av^ 
is an odd function of ni. In the weak- and strong-coupling 
limits we can write down simple expressions for Av^ (see 
appendix B.2): 

At;c«-^l^(l-(An/2)2)3/2 (C/«2i|62) 

Avc ~ —U{1 — |An/2|) sgn(An) {U '> 2t) (63) 




FIG. 10 Plot of Avc for different U and 2t = 1. 


These correspond to the 1st and 4th panels in Fig. 10. 
For small f/, it is of order U‘^ (see appendix B), and has 
little effect. As U increases, it becomes proportional to 
f/, and becomes almost linear in f/, with a large step near 
ni = 1. If we now compare this figure with Fig. 9, we 
see that it is simply the derivative of the previous Ec{ni) 
curve, as stated in Eq. (61). 

The self-consistent KS equations, Eqs. (58) and (59), 
have, in this case, precisely the same form as those of 
restricted HE (or mean-field theory), Eqs. (26) and (37), 
but with whatever additional dependence on rii occurs due 
to Avc{ni). When converged, the ground-state energy is 
found simply from: 

E{ni) = Ts{ni) + Vxt(ni) + Unini) + £'xc(ni)- (64) 

The energy can alternatively be extracted from the KS 
orbital energy via Eq. (16): 

E = 2es + {Ec + AvcAn/2 - Ehx), (65) 

where the second term is the double-counting correction. 
But note the crucial difference here. We consider HE an 
approximate solution to the many-body problem whereas 
DET, with the exact correlation function(al), yields the 
exact energy and on-site occupation, but not the exact 
wavefunction. 


4. THE FUNDAMENTAL GAP 

Now that we have carefully defined what exact KS DFT 
is for this model, we immediately apply this knowledge 
to investigate a thorny subject on the border of many- 
body theory and DFT, namely the fundamental gap of a 
system. 
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A. Background in real space 

Begin with the ionization energy of an A/'-electron sys¬ 
tem: 


I = E{N-1)-E{N), ( 66 ) 

is the energy required to remove one electron entirely from 
a system. We can then define the electron affinity as the 
energy gained by adding an electron to a system, which is 
also equal to the ionization energy of the {N l)-electron 
system: 


A = E{N)-E{NeI). (67) 

In real-space, I and ^4 > 0. For systems which do not 
bind an additional electron, such as the He atom, A = 0. 
The charge, or fundamental, gap of the system is then 


Eg = I-A, ( 68 ) 

and for many materials. Eg can be used to decide if 
they are metals {Eg = 0) or insulators {Eg > 0)(Kohn, 
1964). The spectral function of the single-particle Green’s 
function has a gap equal to Eg. For Coulombic matter. 
Eg has always been found to be non negative, but no 
general proof has been given. 

Now we turn to the KS system of the Welectron system. 
We denote the highest occupied (molecular) orbital as 
eHOMO lowest unoccupied one as Then 

the DFT version of Koopmans’ theorem (Almbladh and 
von Barth, 1985; Almbladh and Pedroza, 1984; Capelle 
et al.^ 2010; Perdew and Levy, 1983; Perdew et a/., 1982; 
Sham and Schhiter, 1983) shows that 

gHOMO ^ 


by matching the decay of the density away from any finite 
system in real space, in the interacting and KS pictures. 
However, this condition applies only to the HOMO, not 
to any other occupied orbitals, or unoccupied ones. In 
particular, the LUMO level is not at —A, in general. 
Define the KS gap as 

E^, = eLUMO_^HOMO^ (70) 

Then Egg does not match the true gap, even with the 
exact XC functional (Baerends et a/., 2013; Sagvolden and 
Perdew, 2008). We write 


Eg — Egs + Axe 

where Axe 7^ 0? and is called the derivative discontinuity 
contribution to the gap (for reasons that will be more 
apparent later) (Perdew, 1985, 1986b). In general. Axe 
appears to always be positive, i.e., the KS gap is smaller 
than the true gap. In semiconductors with especially 
small gaps, such as germanium, approximate KS gaps 
are often zero, making the material a band metal, but 


an insulator in reality. The classic example of a chain of 
H atoms becoming a Mott-Hubbard insulator when the 
bonds are stretched is demonstrated unambiguously in 
Ref. (Stoudenmire et al.^ 2012). 

While this mismatch occurs for all systems, it is es¬ 
pecially problematic for DFT calculations of insulating 
solids. For molecules, one can (and does) calculate the gap 
(called the chemical hardness in molecular systems (Parr 
and Yang, 1989)) by adding and removing electrons. But 
with periodic boundary conditions, there is no simple way 
to do this for solids. Even with the exact functional, the 
KS gap does not match the true gap, and there’s no easy 
way to calculate Eg in a periodic code. In fact, popular 
approximations like LDA and GGA mostly produce good 
approximations to the KS gap, but yield Axe = 0 foi* 
solids. Thus there is no easy way to extract a good approx¬ 
imation to the true gap in such DFT calculations. The 
standard method for producing accurate gaps for solids 
has long been to perform a GW calculation (Aryasetiawan 
and Gunnarsson, 1998), an approximate calculation of the 
Green’s function, and read off its gap. This works very 
well for most weakly correlated materials (van Schilfgaarde 
et al.^ 2006). Such calculations are now done in a variety 
of ways, but usually employ KS orbitals from an approxi¬ 
mate DFT calculation. Recently, hybrid functionals like 
HSE06(Heyd et a/., 2003) have been shown to yield accu¬ 
rate approximate gaps to many systems, but these gaps 
are a mixture of the quasiparticle (i.e., fundamental) gap, 
and the KS gap. Their exchange component produces 
the fundamental gap at the HF level, which is typically a 
significant overestimate, which then compensates for the 
‘too small’ KS gap. While this balance is unlikely to be 
accidental, no general explanation has yet been given. 


B. Hubbard dimer gap 



FIG. 11 Plot of —A, —I, ^lumo ^ function of 

Av with 1/ = 1 and 21 = 1. 
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FIG. 12 Plot of —A, —I, and as a function of 

Av with 1/ = 5 and 21 = 1. 

For our half-filled Hubbard dimer, we can easily cal¬ 
culate both the N ± 1-electron energies, the former via 
particle-hole symmetry from the latter (Carrascal and Fer¬ 
rer, 2012). In Fig. 11, we plot ——^4, and 

eLUMO for [/ = 1 when 2t = 1, as a function of Av. We 
see that A (and even sometimes I) can be negative here. 
(This cannot happen for real-space calculations, as elec¬ 
trons can always escape to infinity, so a bound system 
always has A > 0.) The HOMO level is always at — / 
according to Eq. (69) but the LUMO is not at —A. Here 
it is smaller than —A, and we find this result for all values 
of U and Av. The true gap is / — A, but the KS gap is 
eLUMO ^ which is always smaller. Thus Axe > 0, just 
as for real systems. 

Fig. 11 is typical of weakly correlated systems, where 
Axe is small but noticeable. In Fig. 12, we repeat the 
calculation with 1/ = 5, where now EgEgg at Av = 0, 
but we still see the difference become tiny when Av > U. 
In both figures. Axe is the difference between the red line 
and the green dashed line. In all cases. Axe ^ 0? and this 
has always been found to be true in real-space DFT, but 
has never been proven in general. 


C. Green’s functions 

To end this section, we emphasize the difference between 
the KS and many-body approaches to this problem by 
calculating their spectral functions(Onida et al.^ 2002). 
We define the many-body retarded single-particle Green’s 
function as 

l^-o) (72) 

where label the site indices, a, a' the electron spins, 
and {^4, B} = AB + BA. For the Hubbard dimer at 
N = 1 and 3, |To) is a degenerate Kramers doublet and 
we choose here the spin-j" partner. Fourier transforming 


into frequency, we find for the diagonal component: 


— Gllcrcr(cj) — 

a 


co + E^ - +iS 


+E 




'-E^ + E, 


N-l 



where = (V’aT IV’^)> -^la = ii’a ^\cia\'IPo), 
and (5 > 0 is infinitesimal. Here, a runs over all states of 
the N ± 1-particle systems. The other components have 
analogous expressions. From any component of G, we 
find the corresponding spectral function 


A{uj) = —^G{uj)/7r (74) 


We represent the spectral function (5-function poles with 
lines whose height is proportional to the weights. Via a 
simple sum-rule (Fetter and Walecka, 1971), the sum of all 
weights in the spin-resolved spectral function is 1. There 
are four quasi-particle peaks for V = 2. These peaks are 
reflection-symmetric about uj = Uf2 for the symmetric 
dimer. 

We also need to calculate the KS Green’s function, 
Gs(co’). This is done by simply taking the usual definition, 
Eq. (72), and applying it to the ground-state KS system. 
This means two non-interacting electrons sitting in the 
KS potential. The numerators vanish for all but single 
excitations. Thus the energy differences in the denomi¬ 
nators become simply occupied and unoccupied orbital 
energies. Since there are only two distinct levels (the 
positive and negative combinations of atomic orbitals), 
there are only two peaks, positioned at the HOMO and 
LUMO levels, with weights: 




1 

2 



At>s/2 \ 

^{Avs/2y+ty ’ 


(KS) (75) 


and the sign between the contributions on the right is 
negative in the L term. Thus the symmetric dimer has 
KS weights of 1/2. 

In Fig. 13 we plot the spectral functions for the symmet¬ 
ric case, for U = 1, when 2t = 1. Each pole contributes a 
delta function at a distinct transition frequency, which is 
represented by a line whose height represents the weight. 
The sum of all such weights adds to 1 as it should, and 
the peaks are reffection-symmetric about 17/2 = 0.5. The 
gap is the distance between the highest negative pole (at 
—I) and the lowest positive pole (at —A). We see that 
the MB spectral function also has peaks that correspond 
to higher and lower quasi-particle excitations. If we now 
compare this to the exact KS Green’s function Gg, we 
see that, by construction, Gg always has a peak at — 
whose weight need not match that of the MB function. 
It has only two peaks, the other being at ^lumo^ which 
does not coincide with the position of the MB peak. This 
is so because the KS scheme is defined to reproduce the 
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FIG. 13 Spectral function of symmetric dimer for U — 1, 
Av = 0, and 2t = 1. The physical MB peaks are plotted in 
blue, the KS in red. Here I — 0.1, A = —1.1, and = 0.9, 

corresponding to Av = 0 in Fig. 11. 

ground-state occupations, nothing else. But clearly, when 
U is sufficiently small, it is a rough mimic of the MB 
Green’s function. The larger peaks in the MB spectral 
function each have KS analogs, with roughly the correct 
weights. One of them is even at exactly the right position. 
Thus if a system is weakly correlated, the KS spectral 
function can be a rough guide to the true quasiparticle 
spectrum. 



FIG. 14 Same as Fig. 13, but now 1/ = 5. Here I = —0.3, 
A = —4.7, and = 1.3, corresponding to Av = 0 in Fig. 

12 . 

On the other hand, when U 2 the KS spectral 
function is not even close to the true MB spectral function, 
as illustrated in Fig. 14. Now the two lowest-lying MB 
peaks approach each other, as do the two highest lying 
peaks, therefore increasing the quasi-particle gap. In 
addition, the weights tend to equilibrate with each other. 
In fact, when U ^ oo and/or t —> 0, those two lowest- 
lying peaks gather together at u; = 0, having both the 
same weight of 1/4. And similarly the two highest-lying 


peaks merge at cu = 7/, also with a weight of 1/4. They 
become the precursors of the lower and upper Hubbard 
bands with a quasi-particle gap equal to U. If more sites 
are added to the symmetric dimer, other quasi-particle 
peaks appear, that also merge into the lower and upper 
Hubbard bands as 7/ ^ oo. Notice that the spectral 
function has significant weights for transitions between 
states that differ from the HOMO and LUMO, and are 
forbidden in the KS spectral function for large U. In 
Fig. 14, we see that not only there is a large difference 
between the gaps in the two spectral functions, but also 
the KS weights are not close to the MB weights. The 
only ‘right’ thing about the KS spectrum is the position 
of the HOMO peak. 



FIG. 15 Same as Fig. 13, but now 1/ = 1, Av = 2. Here 
I = 0.27, A = —1.27, and = 1.25, corresponding to 

Ai; = 2 in Fig. 11. 

In Fig. 15, we plot the spectral functions for Av = 2 
and U = 1, to see the effects of asymmetry on the spectral 
function. Now the system appears entirely uncorrelated, 
and the KS spectral function is very close to the true 
one, much more so than in the symmetric case. Here Axe 
is negligible. The asymmetry of the potential strongly 
suppresses correlation effects. In Fig. 16, we see that the 
effects of strong U are largely quenched by a comparable 
Av. Here Axe is small compared to the gap, but not all 
KS peak heights are close to their MB counterparts. 

The situation is interesting even for the ‘simple’ case, 
A" = 1, in which the ground-state is open-shell(Gritsenko 
and Baerends, 2004). Here the interacting spin-t and 
Green’s functions differ. To understand why, we choose 
the N = 1 ground state to have spin t- This state 
has energy ^(1) = + (Au/2)^. Adding a ^-spin 

electron takes the system to the different singlet states at 
A = 2, and to the triplet state with Sz = 0. One of them 
is the ground state at A = 2 whose energy E{2) < 0 is 
given in Eq. (Al) in the appendix. In contrast, adding an 
t-spin electron takes the interacting system to the triplet 
A = 2 state with Sz = whose energy is trivially given 
by E'(2)trip = 0. Annihilating an ^-spin electron takes the 
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FIG. 16 Same as Fig. 15, but now 1/ = 5, Av = 5. Here 
I = —1.8, A = —3.2, and corresponding to Av = 5 

in Fig. 12. 

system to the vacuum, while it is impossible to annihilate a 
|-spin electron. These clearly illustrates that the number 
and energy of the poles in and is different: G^ 
has only two quasi-particle peaks, with trivial energies 
F;(2)trip - ^(1) = and E{1) - E{0) = 

— + (Au/2)^. This last expression corresponds to 

the ionization energy I = E{0) — E{1) = + (Au/2)^. 

Gi has four quasiparticle peaks, all corresponding to 
adding a ^-spin electron, with non-trivial energies. The 
lowest of these corresponds to the electron affinity A = 
E{1) — E{2) = — E{2). In other words, 

ionization involves either removing an ^-spin electron 
(hence seen as a pole in G^) or adding a ^-spin electron 
(hence seen as a pole in G^). The interacting gap is 
Eg = I-A = 2 + (Av/2)2 + E{2). 

We turn now to the KS Green’s function. For N = 
the KS on-site potentials equal the true on-site potentials, 
±Av/2. So the ground-state (chosen again to have spin 
t) has energy ^s(l) = —+ (Aus/2)^. Since the other 
state has energy ^s(l): a-^d a second ^-electron occupies 
that state, the total KS energy is F^(2)s^=i = 0. On the 
other hand, annihilating the t electron costs an energy 
E{1). This shows that the t-spin KS and interacting 
Green’s functions are identical to one other and trivial 
for AT = 1. Thus / = -gHOMO ^ ,yW+JKv/^. This 
result is specific to this model. 

Removing a ^-spin KS electron is impossible, just 
as in the interacting case. However, adding it means 
having either two opposite-spin KS electrons with 
the same energy + (Aus/2)^, or having one 

with energy —\/E (Aus/2)^ and another with energy 

^/E + (Aus/2)^. The first case corresponds to the KS 
ground-state with energy —2 y^E + (Aus/2)^, while the 
second one is an excited state with energy 0. The KS 
value for the electron affinity is As = ^s(l) ~ £'s(2) = 
y/E + (Aus/2)^, which differs from the interacting value. 


Furthermore, the KS gap Egs = 0 is clearly an incorrect 
estimate of the true interacting gap, which is given by 



CO 


FIG. 17 Spin-| resolved spectral function for iV = 1 and 
U = 1, Av = 2. Here I = 1.12, A = 0.27, and = 

eHOMO ^ ^2. 



FIG. 18 Spin-j. resolved spectral function for N = 1 and 

= 5, Aa = 2. Here I = 1.12, A = -0.90, and = 

eHOMO ^ 12. 

Figs. 17 and 18 show the spectral function associated 
with Gi for the many-body and KS Green’s functions for 
N = 1 and Av = 2. In the first, U = 1, so it is relatively 
asymmetric, whereas in the second, C/ = 5, making it 
close to symmetric. Thus the HOMO is at the lowest 
red line, and matches exactly the LUMO, with a KS gap 
of zero. Thus Axe the gap of the interacting system. 
We see that in the first figure, correlations are weak and 
the KS spectral function mimics the physical one, but in 
the second figure {U = 5), they differ substantially, even 
though = 1! 

The difference in expressions for spin species is illus¬ 
trated further by work analyzing Koopmans’ and Janak’s 
theorems for open-shell systems (Gritsenko and Baerends, 
2002, 2004; Gritsenko et a/., 2003; Griining et a/., 2003). 
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Self-energy approximations beyond GW have been per¬ 
formed on the Hubbard dimer (Romaniello et a/., 2012 , 
2009), as well as a battery of many-body perturbation 
theory methods (Olsen and Thygesen, 2014) though only 
for the symmetric case. 

The bottom line message of this subsection is that the 
KS spectral function does not match the quasiparticle 
spectral function, because it is not supposed to. However, 
the main features of a weakly correlated system are loosely 
approximated by those of the KS function, with the gap 
error shifting the upper part of the spectrum relative 
to the lower part. This is the motivation behind the 
infamous scissors operator in solid-state physics. A very 
accurate DFT approximation can (at best) approximate 
the KS spectral function, not the many-body one. The 
exact XC functional does not reproduce the quasiparticle 
gap of the system. For strongly correlated systems, there 
are often substantial qualitative differences between the 
MB and KS spectral functions. These are some of the 
limitations of KS-DFT. that, e.g., DMFT is designed to 
overcome (Georges et al.^ 1996). 

5. CORRELATION 


In Figs. 8 and 19, we plot both and Tc, respectively, for 
several values of U (with 2t = 1). When U is small, Tc ~ 
—Tc. However, for U ^ 2t, we see that although 
becomes very large (in magnitude), Tc remains finite and 
in fact, Tc never exceeds 2 t as proven in Appendix C. We 
can define a measure of the nature of the correlation (Burke 
et a/., 1997): 

Pcovv = 1^ I 5 (80) 

As U —y 0, Pcorr —^ I 5 while as U —y CXD, Pcorr —^ 0. 
Thus /3corr close to 1 indicates weak correlation, /3 small 
indicates strong correlation. We plot Peon as a function 
of U for several values of Av in Fig. 20. Although Peon 
is monotonically decreasing with U for Av = 0 , we see 
that the issue is much more complicated once we include 
asymmetry. The curve for each Av remains monotonically 
decreasing with U. But consider U = 2 and different 
values of Av. Then Peon at first decreases with Au, i.e. 
becoming more strongly correlated, but then increases 
again for Av > 1/, ultimately appearing less correlated 
than Av = 0 . 


A. Classifying correlation: Strong, weak, dynamic, static, 
kinetic, and potential 

There are as many different ways to distinguish weak 
from strong correlation as there are communities that 
study electronic structure. Due to the limited degrees of 
freedom (namely, one), these all overlap in the Hubbard 
dimer. We will discuss each. 

The most important thing to realize is that correlation 
energy comes in two distinct contributions: kinetic and 
potential. These are entirely well-defined quantities within 
KS-DFT. The kinetic correlation energy is: 

Tc = T - Ts (76) 

for a given density. Note that we could as easily call this 
the correlation contribution to the kinetic energy. The 
potential correlation energy is: 

l/c = Kee-^HX, (77) 




FIG. 19 Plot of exact Tc (blue line) and Tc,par (red dashed 
line) for different U and 21 = 1. 


and could also be called the correlation contribution to 
potential energy. For future notational convenience, we 
also define = Tx, i.e., there is no kinetic contribution 
to exchange. Then, from Eq. (18), we see 

Tc=Tc + t/c. (78) 

We can now use these to discuss the differences be¬ 
tween weak and strong correlation. First note that, by 
construction, and as shown for our dimer in appendix C, 

(79) 


Quantum chemists often refer to dynamic versus static 
correlation. Our precise prescription in KS-DFT loosely 
corresponds to their definition, replacing dynamic by 
kinetic, and static by potential. Thus, considering an 
H 2 molecule with a stretched bond, the Hubbard model 
applies. As the bond stretches, t vanishes, and U/2t grows. 
Thus Peon ^ 0 as T ^ 00 . The exact wavefunction, the 
Heitler-London wavefunction (Heitler and London, 1927), 
has only static correlation in this limit. In many-body 
language, it is strongly correlated. In DFT language, the 
fraction of correlation energy that is kinetic is vanishing. 


Tc <0, Tc > 0, Uc < 0. 
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FIG. 20 Plot of y^corr = Tc/\Ec\ as a function of U with 
2t = 1. 


B. Adiabatic connection 

With the various contributions to correlation well- 
defined, we construct the adiabatic connection (AC) for¬ 
mula (Gunnarsson and Lundqvist, 1976; Langreth and 
Perdew, 1975) for the Hubbard dimer. The adiabatic 
connection has had enormous impact on the field of DFT 
as it allows both construction (Adamo and Barone, 1999; 
Becke, 1993a,b; Ernzerhof and Scuseria, 1999; Perdew 
et al.^ 1996b; P.Mori-Sanchez et a/., 2006), and under¬ 
standing (Burke et a/., 1997; Peach et a/., 2008; Perdew 
et a/., 1996b), of exact and approximate functionals solely 
from their potential contributions. 

In many-body theory, one often introduces a coupling- 
constant in front of the interaction. In KS-DFT, a cou¬ 
pling constant A is introduced in front of the electron- 
electron repulsion but, contrary to traditional many-body 
approaches, the density is held fixed as A is varied (usually 
from 0 to 1). Via the Hohenberg-Kohn theorem, as long 
as there is more than 1 electron, this implies that the 
one-body potential must vary with A, becoming 
By virtue of the density being held fixed, i;^^^(r) = 'i;s(r) 
while = v{r). Thus A interpolates between the 

KS system and the true many-body system. Addition¬ 
ally, A ^ oc results in the strictly correlated electron 
limit (Gori-Giorgi and Seidl, 2010; Liu and Burke, 2009a; 
M. Seidl and Levy, 1999; Malet and Gori-Giorgi, 2012; 
Seidl et a/., 2007) which provides useful information about 
real systems that are strongly correlated. 

The adiabatic connection for the Hubbard dimer is very 
simple. Define the XC energy at coupling constant A by 
simply multiplying U hy X while keeping An fixed: 

([/, An) = E^aiXU, An). (81) 

Application of the Hellman-Feynman theorem(Feynman, 
1939) yields (Gunnarsson and Lundqvist, 1976; Harris and 


Jones, 1974; Langreth and Perdew, 1975, 1977): 

dE^^{XU,An) Vxc(AV,An) 

- Tx ---A-’ 

where An) is the potential contribution to the XC 

energy, i.e., = E^ and 


Uc{XU) = Vee(AV) -XE^^iU). (83) 


Thus, we can extract Tc solely from our knowledge of 
Ec{U) via 


Tc = Ec - Vc = Ec - ^ . (84) 

Thus, any formula for Ec, be it exact or approximate, 
yields a corresponding result for Tc and Vc, and vice 
versa(D. Frydel and Burke, 2000). We may then write 

E^c{U,An) = — C/xc(A?7,An), (85) 

and this is the infamous adiabatic connection formula 
of DFT(Gunnarsson and Lundqvist, 1976; Langreth and 
Perdew, 1975). We denote the integrand as Vc(A), defined 
as 


„ l7„{Ai7) dE^(XU) 

C/c(A) = (86) 

Plots of Uc{X) from Eq. (86) are called adiabatic connec¬ 
tion plots, and can be used to better understand both 
approximate and exact functionals. In Fig. 21, we plot a 
typical case for 1/ = 2t and Av = 0. They have the nice 
interpretation that the value at A = 1 is the potential 
correlation energy, Vc, the area under the curve is Ec, 
and the area between the curve and the horizontal line at 
f/c(l) is — Tc- Furthermore, one can also show(Levy and 
Perdew, 1985) 


dU^cW 

dX 


(87) 


from known inequalities for Tc(A) and Ec{X). This is 
proven for our problem in appendix C. Interestingly, such 
curves have always been found to be convex when ex¬ 
tracted numerically for various systems (Fuchs et a/., 2005; 
Puzder et a/., 2001), but no general proof of this is known. 
The Hubbard dimer also exhibits this behavior. A proof 
for the dimer might suggest a proof for real-space DFT. 

In Fig. 21 we plot Uc{X)/U for Av = 0 and Av = 2, 
with various values of U. From the above formulas, one 
can deduce that the area between the curve and the 
horizontal line at Vc(l) is —Tq. Thus as U grows, the 
curve moves from being almost linear to decaying very 
rapidly, and (3 varies from 1 down to 0. 

In Fig. 21, we show U up to 10 (for 2t = 1), to show the 
effect of stronger correlation. Not only has the magnitude 
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of the correlation become larger, but the curve drops more 
rapidly toward its value at large A. pcow — 0-9 for Ai; = 0 
and 1/ = 1, but /3corr — 0.2 for /\v = 0 and U = 10, 
reflecting the fact that the increase in correlation is of the 
static kind. 

The weakly correlated limit has been much studied 
in DFT. Perturbation theory in the coupling constant 
is called Goerling-Levy perturbation theory (Gorling and 
Levy, 1993). For small A, 

Ua{XU) = + X^U^c^ + ... (A ^ 0). (88) 

In Appendix B.2, we show that 

and 

for the dimer. This yields, for Tc, 

To = -\x^U^c^ - 

showing that /3 ^ 1 as 1/ (or A) vanishes. For any system, 
determines the initial slope of Uc{X). 

On the other hand, in the strongly correlated limit, in 
real-space(Gori-Giorgi et a/., 2009; Liu and Burke, 2009a). 

Eq A(5o + A + A ^^ 2 ...), (A oo) (92) 

where (/c = 0,l,2...) are coupling-invariant functionals 
of n(r)(Liu and Burke, 2009b). The dominant term is 
linear in U. Physically, it must exactly cancel the Hartree 
plus exchange contributions, since there is no electron- 
electron repulsion to this order when each electron is 
localized to separate sites. Correctly, such a term cancels 
out of Tc, so that its dominant contribution is 0(1). 
From Appendix B.2, we see that the Hubbard dimer has 
a different form, involving only integer powers of A: 

E(j A Bq + Bi + B2 / A-f-... (A^oc) (flfl) 

where 

Bo(An) = -U{1 - Anl2fl2, (94) 


Bi{An) = 2t^/l - An/2 (^/l + An/2 - \/An), (95) 

and 

B 2 {An) = {1 - An/2)t^/U. (96) 


But both this term and the next cancel in the total energy 
(at half Ailing), so that the ground-state energy is 0(1/t/), 
i.e., extremely small as U grows: 




This illustrates that, although the KS description is exact, 
it becomes quite contorted in the large U limit. This 
has been implicated in convergence difficulties of the KS 
equations, even with the exact XC functional, because 
the KS system behaves so differently from the physical 
system(Wagner et al.^ 2014). 



FIG. 21 Adiabatic connection integrand divided by U for 
various values of U. The solid lines are Av = 2 and the dashed 
lines Av = 0. Asymmetry reduces the correlation energy but 
increases the fraction of kinetic correlation. 


6. ACCURATE PARAMETRIZATION OF CORRELATION 
ENERGY 

Although the Hubbard dimer has an exact analytic 
solution when constructed from many-body theory, the 
dependence of F{An) (or equivalently Tc(An)) is only 
given implicitly. While this is technically straightforward 
to deal with, in practice it would be much simpler to 
use if an explicit formula is available. In this section, we 
show how the standard machinery of DFT can be applied 
to develop an extremely accurate parametrization of the 
correlation energy functional. 

An arbitrary antisymmetric wavefunction is character¬ 
ized by 3 real numbers where |12) means an electron at 
site 1 and site 2, etc.: 

1^) = a (|12) + |21)) + /3i 111) + |22). (98) 

Normalization requires 20^ + terms of 

these parameters, the individual components of the energy 
are rather simple: 

T = -4ta(A+^2) 

V,, = U{pl+pl) 

V = -Av{pl -Pl), (99) 

SO that the variational principle may be written as 

E= min E{a,j3i,p2)‘ (100) 

a,/3i,/32 

l=2cx^+(3l+(3l 


( 97 ) 
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The specific values of these parameters for the ground- 
state wavefunction are reported in appendix A. 

For this simple problem, we are fortunate that we can 
apply the Levy-Lieb constrained search method explicitly. 
A variation of this method was used for the derivation of 
the exact functional of the single- and double-site Ander¬ 
son model and the symmetric Hubbard dimer (Carrascal 
and Ferrer, 2012). The functional F[n] is defined by 
minimizing the expectation value of T -t- over all pos¬ 
sible wavefunctions yielding a given n(r). In real-space 
DFT, there are no easy ways of generating interacting 
wavefunctions for a given density. But here. 

An = 2{pl - pl), (101) 


These forms are chosen so qq is exact to second- and first- 
order in the weak- and strong-coupling limits respectively, 
and to first- and second- order in the symmetric and 
asymmetric limits respectively. Use of this to construct 
an approximation to F, f{go{p)^p), yields very accurate 
energetics. The maximum energy error, divided by U, is 
0 . 002 . 

But for some of the purposes in this paper, such as 
calculations of Tc, even this level of error is unacceptable. 
We now improve on go{p) using the adiabatic connection 
formula of Sec 5.B. Like F, we can define functions of two 
variables for each of the correlation components. Write 

ec(fl, p) = fig, P) - Tsip) - -E'hx(p)- (109) 


which allows us to simply eliminate a parameter, e.g., p 2 
in favor of An. Thus 

F[An] = min [Tia, pi, An)+ Veeia, pi, An)]. 

a2+/3j = i(l+l^) 

( 102 ) 

With normalization and the density constraint, only one 
parameter is left free. There exist several possible choices 
for this. If we choose g = 2a {/3i + ^ 2 ) which corresponds 
to the hopping term, then after some algebra the func- 
tion(al) can be written nicely as 

F{p) = min f{p,g) (103) 

9 

with the intermediate quantity 

fip,g) =-“^tg+ Uhig,p), (104) 


where Tg and Fhx are from Eqs. (54) and (53), respec¬ 
tively. The kinetic and the potential correlation are given 
by 

tGig,p) = T -Ts = -2t (g- i/l- (110) 

Ucig,p) = Vee - £^hx = u [hig,p) - (1 + y)/2](;Lll) 

and their sum yields edg^p)- If we insert g{p), the exact 
minimizer of /(^, p), into any of these expressions, we get 
the exact answers. 

But recall also that one can extract Uc from the deriva¬ 
tive of Fc with respect to the coupling constant A, i.e., 

Uc = dEdX)/dX\x=i. ( 112 ) 

Now for any g and ec(p), we can find the A dependence 
by replacing U by XU. Thus 


and 


Kg, p) 


g^ (1 - \/i - - p^) + 2/9^ 

2(^2 + p2) 


(105) 


Note that both t and U appear linearly in f{g^p). The 
minimization yields a sextic polynomial, equation (Bl), 
that g must satisfy. The weak-coupling, strong-coupling, 
symmetric, and asymmetric limits of g are given in ap¬ 
pendix B. 

Our construction begins with a simple approximation 
to gip): 


d-ecig,A) _ dE^iX) ^ dE^iX) dg 

d\ “ d\ dg ^ ^ ^ 

Since Tg and Fhx do not depend on p, the minimization 
of / reduces to de^jdg = 0, so for the exact g the second 
term on the right of Eq. (113) is always zero. But it does 
not vanish for po- 

Equating Eqs. (112) and (113) and using the definitions, 
we find the following self-consistent equation for g: 


^ Id^ dg{X) 
2t 2t dg dX 


(114) 


goip) 

where 


> (1 - p) (1 + p (1 + (1 + p)^uai{p, u))) 

1 + (1 + p)^ua2ip,u) 


(106) 


aiip,u) = anip) +u ai 2 (p), 


(107) 


We may use this to improve our estimate for g. Simply 
evaluate the right-hand side at po, to find: 


dh dg{X) 


A=1 


( 115 ) 


where 


and 


a2i = ^\/(l - p)p/2, an = a2i(l+ P ^), 
ai2 = 


dgjX) (l-p)(l + p)3 

9 ^ so 2po(l + (1 + p)^a2(A))2 

^ [/^(l + (1 + /5)^®2(A))ai(A) 
— (1 + p(l + p)^ai(A))a2(A)]. 


022 = 012 / 2 . 


(108) 
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FIG. 22 Error in Ec,par{p)/U for different U and 2t = 1. 


The new Fpar and £’c,par are then obtained by using gi 
in Eqs. (104) and (109). Using gi, dEc^pav/dg 7 ^ 0 still, 
but the error with gi is much lower than with g^. We 
plot the the relative error, {E^ — Ec,par)/U for several U 
in Fig. 22. The maximum relative error is reduced by 
almost two orders of magnitude (from 2 x 10 “^ to 5 x 10 “^) 
in the region U ~ 2 — 6 , An ~ 0.25, where go has the 
largest error. The other regions are also improved. For 
(Tc — T'c,par)/U and (Uc — Uc,par)/U the improvement is 
just of one order of magnitude (from 2 x 10 “^ to 2 x 10 “^ 
in both cases relative to the maximum), with different 
sign, so there is an error cancellation that yields the larger 
reduction of the E^ error. We anticipate that g could be 
improved even further by iteration. 


in the KS scheme to calculate the correlation energy of 
our Hubbard dimer self-consistently. If our parametriza- 
tion were perfect, we would recover the exact densities 
and energies from our KS calculation without having to 
solve the many-body problem. These are plotted in Figs. 
23, together with the absolute errors committed by the 
parametric function(al). Notice that in Figs. 8 and 19 
the results obtained from the parametric function(al) are 
indistinguishable from the exact results. We recommend 
the use of go for routine use, and gi for improved accuracy. 
We hope the methodology developed here might prove 
useful to improve accuracy of correlation functionals in 
other contexts(Site, 2015). 

We can define the starting point of our parameterization 
in a multitude of ways. In this section we defined it such 
that the parameter corresponds to the hopping term. 
Another possible choice favors the electron-electron term. 
Define, 

f 2 {f,p) = - 2^71 -/ (^Vf + p+Vf-p) + uf. 

(116) 

Another choice captures the asymmetric limit. Define, 

72 2 

= (117) 

Then, 

Fi.P) = min/2(/, p) = rainfsH, p) (118) 

J ^ 

These also yield high order polynomial equations when 
minimized. The present parameterization, Eq. (106), is 
quantitatively superior for nearly all values of U, and Av 
of interest. 

7. APPROXIMATIONS 



xio'' U=5 




FIG. 23 Top row: Error in density as a function of Av. 
Bottom row: Error in ground-state energy as a function of Av 
and 21 = 1. 


To test the validity of our parametrization, we use it 


The usefulness of KS-DFT derives from the use of ap¬ 
proximations for the XC functional, not from the exact 
XC which is usually as expensive to calculate as direct 
solution of the many-body problem (or more so). While 
the field of real-space DET is deluged by hundreds of dif¬ 
ferent approximations (Marques et al.^ 2012a), (relatively 
few of which are used in routine calculations (Pribram- 
Jones et a/., 2015)) few approximations exist that apply 
directly to the Hubbard dimer. The two we explore here 
are illustrative of many general principles. 


A. Mean-field theory: Broken symmetry 

Since time immemorial, or at least the 1930’s, folks 
have realized the limitations of restricted HF solutions 
for strongly correlated multi-center problems, and per¬ 
formed broken-symmetry calculations (Coulson and Fis¬ 
cher, 1949). For example, in many-body theory, Ander¬ 
son solved the Anderson impurity model for a magnetic 
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atom in a metal (Anderson, 1961) by allowing symmetry 
breaking, several years before Kondo’s ground-breaking 
work(Kondo, 1964). In quantum chemistry, Coulson 
and Fischer identified the Coulson-Fischer point of the 
stretched H2 molecule where the broken symmetry solu¬ 
tion has lower energy than the restricted solution(Coulson 
and Fischer, 1949). Modern quantum chemists like to spin 
purify their wavefunctions, but DFT hardliners (Perdew 
et al.^ 1995) claim the broken-symmetry solution is the 
‘correct’ one (for an approximate functional). The exact 
KS functional, as shown in all previous sections, yields 
the exact energy and spin densities, while remaining in a 
spin singlet. 

If we do not impose spin symmetry, the effective poten¬ 
tial in mean-field theory becomes (Sec 2.B): 

=Vi + U Tlia, (119) 

with cr = +1 for spin up, a = —1 for spin down and 
a = because the change in the effective field is caused 
by the other electron. Writing rii = rii^^ + rrii = 

~ '^2, and defining 


Av^ = Av — ^ (An — a Am), 


U 


( 120 ) 


and 


tf=t^l + {Avf/2ty, (121) 

we find the eigenvalues are: 

TT /eff 

ef!^,=j{N-aM)±^, ( 122 ) 

where A = 2 is the number of particles and M is the 
total magnetization. We find the ferromagnetic solution 
(M = 2) to be everywhere above the antiferromagnetic 
solution (M = 0), and for M = 0: 


E=^il 
2 ^ 


An^ + ( 123 ) 


where Am = 0 is the paramagnetic (spin singlet) solution, 
and corresponds to our original mean-field or restricted 
Hartree-Fock solution. We minimize this energy with 
respect to An and Am, given by 


An = ^ 


Avf 

1^’ 


Am = a 




(124) 


These antiferromagnetic (AFM) self-consistency equations 
always have the trivial solution Am = 0, which corre¬ 
sponds to the restricted MF solution(RHF). However, 
there exists a non-trivial solution Am 7^ 0 for sufficiently 
large values of U. 

In Fig. 24, we plot An and Am for both restricted 
and unrestricted HF solutions for U = b. The solutions 
coincide for large Av^ but below a critical value of Av^ 



Av 


FIG. 24 Plots of An for HF and BALD A as a function of Av 
for U — b and 2t = 1. The crossover from the charge-transfer 
to the Mott-Hubbard regime happens at F Av. 



FIG. 25 Ground-state energy of the unrestricted Hartree-Fock 
(thick dashed line), restricted Hartree-Fock (dot dashed line), 
and exact ground-state (thin solid line) of the Hubbard dimer 
as a function of Av for several values of U and 2t — 1. The 
dot shows the Goulson-Fischer point at which the symmetry 
breaks spontaneously. For smaller Av the UHF energy is below 
RHF while for larger Av they are the same. 


they differ. The UHF solution has a significantly lower 
An, which is much closer to the exact An. 

In Fig. 25, we plot the energies, showing that the UHF 
solution does not rise above zero, and mimics the exact 
solution rather closely. For large F, at ni = 1, we can 
compare results analytically: 

U 2t^ 4t^ 

E^--2t (RHF), (UHF), (exact) 

(125) 

confirming that the UHF energy is far more accurate 
than the RHF energy, and recovers the dominant term in 
the strongly correlated limit. Note that the symmetric 
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case is atypical: The constant terms vanish, both exactly 
and in UHF, so the leading terms is O(!/[/), and its 
coefficient in UHF is underestimated by a factor of 2. 
The slope of the exact result is two times larger than 
UHF. Of course, the exact solution is a spin-singlet, so 
the symmetry of the UHF solution is incorrect, but its 
energy is far better than that of RHF. This is called the 
symmetry dilemma in DFT(Perdew et a/., 1995): Should 
I impose the right symmetry at the cost of a poor energy? 
Note that the exact KS wavefunction is also a singlet, 
so a broken-symmetry DFT solution produces the wrong 
symmetry for the KS wavefunction. 


B. BALDA 


Verdozzi, 2008; Verdozzi et a/., 2011) and also transport 
properties(Kurth et a/., 2010; Vettchinkina et a/., 2013), 
as well as using BALDA as a gateway to calculate time- 
dependent effects in 3-D(Karlsson et a/., 2011a). There 
has been significant interest in using BALDA to under¬ 
stand the derivative discontinuity in both DFT and TD- 
DFT(Kurth et a/., 2010; Xianlong et a/., 2012, 2006c; Ying 
et a/., 2014). Additionally, the BALDA approach has been 
developed for other BA-solvable fermionic lattice systems 
aside from the Hubbard model(Abedinpour et al.^ 2007b; 
Mirjani and Thijssen, 2011; Schenk et a/., 2008; Xianlong 
et a/., 2006a), such as the Anderson model (Bergfield et a/., 
2012; Kurth and Stefanucci, 2013; Liu et a/., 2012), as 
well as bosonic systems (Hao and Chen, 2009; Wang et a/., 
2012; Wang and Zhang, 2013). 


In real-space DFT, the local density approximation 
(LDA) was first suggested by Kohn and Sham(Kohn and 
Sham, 1965), in which the XC energy is approximated 
at each point in a system by that of a uniform gas with 
the density at that point. Another way to think of this is 
that one decides to make a local approximation, and then 
chooses the uniform gas XC energy density to ensure exact¬ 
ness in the uniform limit. On the lattice, we must switch 
our reference system to incorporate Luttinger-liquid corre¬ 
lations instead of Fermi-liquid correlations (Haldane, 1981). 
The infinite homogeneous Hubbard chain plays the role 
of the uniform gas. This can be solved exactly via Bethe 
ansatz(Lieb and Wu, 1968), and the corresponding LDA 
was first constructed and tested in Ref. (Schonhammer 
et a/., 1995). Later, Capelle and collaborators (Capelle 
et al.^ 2003; Franca et a/., 2012b; Lima et al.^ 2002, 2003; 
Xianlong et al.^ 2006c) used the exact Bethe ansatz solu¬ 
tion to create an explicit parameterization for the energy 
per site, and called this Bethe Ansatz LDA, or BALDA. 

Since its inception, BALDA has been applied to many 
different problems including disorder and critical be¬ 
havior in optical lattices(Campo Jr. et a/., 2009; Xian¬ 
long, 2008), spin-charge separation(Vieira, 2012, 2014) 
and effects of spatial inhomogeneity (Lima et al.^ 2007; 
Silva et a/., 2005) in strongly correlated systems, con¬ 
fined fermions both with attractive and repulsive inter- 
actions(Campo and Capelle, 2005), current DFT on a 
lattice (Akande and Sanvito, 2012), electric fields and 
strong correlation (Akande and Sanvito, 2010), and vari¬ 
ous critical phenomena in 1-D systems (Abedinpour et a/., 
2007a; Franca et al.^ 2012a). Extensions to include spin- 
dependence (BALSDA) have been principally used for 
studying density oscillations (Vieira et a/., 2008; Xianlong, 
2012), and fermions in confinement (Hu et al.^ 2010; Xi¬ 
anlong, 2013; Xianlong and Asgari, 2008). A thermal 
DFT approximation on the lattice has been constructed 
using BALDA(Xianlong et a/., 2012). BALDA has also 
been used as an adiabatic approximation in TD-DFT to 
calculate excitations (Khosravi et al.^ 2012; Kurth and 
Stefanucci, 2011; Li et al.^ 2008; Uimonen et al.^ 2011; 


LU 





FIG. 26 Ground-state energy versus Av for several U, with 
2t = 1. The BALDA energies are evaluated self-consistently. 

We use here the semi-analytical approach to 
BALDA(Lima et al.^ 2003; Xianlong et al.^ 2006c) where 
the expressions are given in Appendix D. In Fig. 26 we 
plot the BALDA ground-state energy as a function of 
Av for several values of U. At first glance, it seems to 
do a good job in all regimes. In particular, for either 
very weak correlation {U = 0.2) or very strong correlation 
{U = 100), it is indistinguishable from the exact curves. 
However, for moderate correlation (1 < U ^5) where 
Av < U, it appears to significantly underestimate the 
magnitude of E. 

Even for the strong correlation regime, its behavior is 
not quite correct. For the symmetric case: 

>0 (f7»2t) (126) 

Thus, for Av = 0 and U = 100 in Fig. 26, BALDA is 
in serious error, but this cannot be seen on the scale of 
the figure. The origin of this error is easy to understand. 
BALDA’s reference system is an infinite homogeneous 
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chain, and we are applying it to a finite inhomogeneous 
dimer. The error is in the correlation kinetic energy, 
which comes from the difference between the exact and 
KS kinetic energies. The tight-binding energy for an 
infinite homogeneous chain is different from that of the 
dimer, and this difference is showing up (incorrectly) in 
the correlation energy. 


C. BALDA versus HF 





Av Av 

FIG. 27 Plots of the RMF, UMF, and BALDA API = 
^approx _ ^exact ^ function of Av for U = 0.2, 1, 5, and 10. 
For small U the RMF and UMF results are indistinguishable. 
Here 21 = 1. 

Lastly we compare BALDA and both the restricted 
and unrestricted Hartree-Fock approximations. In Fig. 
27, we plot the errors made in the ground-state energy 
of all three approximations. For U < 1, HF does not 
break symmetry, and so UHF=RHF. For very small U, 
the energy error is comparable to HF. For 1/ = 1, BALDA 
is better than HF. For larger U, UHF produces a lower 
energy than HF, and almost everywhere is more accurate 
than BALDA. The sole exception is at precisely U ~ Au, 
where BALDA is much better. In Fig. 24, we compare 
BALDA and UHF densities to the exact density for 1/ = 5 
as a function of Av. Although BALDA does not have a 
symmetry-breaking point, it unfortunately has a critical 
value of Av where An vanishes incorrectly. This is the 
origin of the cusp-like features in the BALDA energies 
of Figs. 26 and 27. In fact, the BALDA density appears 
somewhat worse than UHF for most Av. But keep in mind 
that the main purpose of BALDA is to produce accurate 
energies without the artificial spin-symmetry breaking of 
UHF. 


8. FRACTIONAL PARTICLE NUMBER 

We will now show a way that one can extract the phys¬ 
ical gap from ground-state DFT. This is done simply by 
changing the number of electrons, but now continuously, 
rather than just at integers. In fact, we already used 
this technology implicitly in Sec 4, but here we make this 
much more explicit. 

A. Derivative discontinuity 

An extremely important concept in DFT is that of the 
derivative discontinuity (Cohen et a/., 2008a; Kurth et a7, 
2010; Mori-Sanchez et al., 2008; Mori-Sanchez and Cohen, 
2014; Mori-Sanchez et al.^ 2009; Perdew and Levy, 1983; 
Perdew et a7, 1982; Schonhammer and Gunnarsson, 1987; 
Yang et al.^ 2012). This is most famous for its implica¬ 
tion for the Kohn-Sham gap of a solid, ensuring that the 
gap (in general) does not match the true fundamental 
(or charge) gap of the solid, as we saw in Sec. 4. The 
expression itself refers to a plot of ground-state energy 
versus particle number N at zero temperature. In seminal 
work(Perdew, 1985; Perdew and Levy, 1983; Perdew et a7, 
1982), it was shown that E{Af) consists of straight-line 
segments between integer values, where A/' is a real vari¬ 
able, where all quantities are now expectation values in a 
grand-canonical ensemble at zero temperature: 

E{Af) = {l-w) E{N) + w E{N + 1), (127) 

and 

«Ar(r) = (1 - w)njv(r)+«;njv+i(r), (128) 

where Af = N + w, i.e., both energy and ground-state 
density are piecewise linear, with a sudden change at 
integer values. 

Then the chemical potential is 

ia = dEtdJ\r=-I {AT <N) 

= -A {Ar>N). (129) 

When we evaluated everything at Y = 2 in Sec. 4, we 
really meant N = 2~. Then Janak’s theorem(Janak, 
1978) shows that, for the KS system, 

IJ. = dE/dAf = {M < N) 

= (V > N) (130) 

This is the proof of the equivalence of I and —gHOMO 
Because the energy is in straight-line segments, the 
slope of E{J\f)^ the chemical potential, jumps dis- 

continuously at integer values. Hence the name, derivative 
discontinuity. The jump in fi across an integer N is then 
Eg = I — the fundamental gap. In the KS system, 
since the energy is given in terms of orbitals and their 
occupations, that jump is simply the KS HOMO-LUMO 

























25 


gap, Egs. Since the KS electrons have the non-interacting 
kinetic energy, and the external and Hartree potentials 
are continuous functionals of the density, the difference 
is an XC effect. Moreover, it implies that jumps by 
this amount as one passes through N^ an integer. 

For solids, addition or removal of a single electron has 
an infinitesimal effect on the density, but the XC discon¬ 
tinuity shifts the conduction band upward by Axe when 
an electron is added, contributing to the true gap. Since 
local and semilocal approximations to XC are usually 
smooth functionals of the density, they produce no such 
shift. They do yield accurate approximations to the KS 
gap of a solid, but not to the gap calculated by adding 
and removing an electron, because of this missing shift. 
Thus we have no general procedure for extracting accurate 
gaps using LDA and GGA. An important quality factor 
in more sophisticated approximations is whether or not 
they have a discontinuity. Orbital-dependent functionals, 
such as exact exchange (EXX in OEP)(Gorling, 2005; 
Krieger et a/., 1992; Kiimmel and Kronik, 2008; Sharp 
and Horton, 1953; Talman and Shadwick, 1976; Yang and 
Wu, 2002) or self-interaction corrected LDA (SIC) (Heaton 
et a/., 1983; Johnson et al.^ 1994; Pederson et al.^ 2014; 
Pemmaraju et a/., 2007; Perdew and Zunger, 1981), often 
capture effects due to the discontinuity quite accurately. 


B. Hubbard dimer near integer particle numbers 



N 


FIG. 28 Plot of E{N) for E = 1, Aa = 0 and 21 = 1. 

In Fig. 28, we plot E{M) for our Hubbard dimer. 
Real-space curves have always been found to be convex, 
although this has never been proven to be generally true. 
The vital part for us is that this equivalence of the HOMO 
level and —I links the overall position of the KS levels to 
those of the many-body system. For fixed particle number, 
only the KS on-site energy difference is determined by 
the need to reproduce the exact site occupancies. But 
this condition also fixes the mean value of the KS on-site 


energy, which in general is non-zero, even though we 
chose the actual mean on-site energy to be zero always. 
In Fig. 2, this is visible in the mean position of the two 
KS on-site potentials. 



FIG. 29 Same as Fig. 11 except with N — 2'^ instead of 
N = 2“. 


Another way to think about this is that function(al) 
derivatives at fixed M leave an undetermined constant in 
the potential, whereas that constant is determined if the 
particle number is allowed to change. We can write many 
equivalent formulas for the discontinuity^: 


Axe — 


5Ex 


dN 


5Ex 


N+ 


N- 


dN 

= v^c{N+)-v^c{N-), 
= Vs{N+)-Vs{N-), 

= es(iV+)-es(iV-), 


(131) 


all of which are true. Thus another way to find the gap 
from a KS system is to occupy it with an extra infinites¬ 
imal of an electron, and note the jump in potentials or 
eigenvalues. To illustrate this, in Fig. 29 we replot Fig. 
11, but now for A = 2+, showing that now the LUMO 
matches —A, and the difference between the HOMO and 
I is Axe* 

In Fig. 30 we plot Axe for N = 2 for various 1/, as a 
function of A'z;, scaling each variable by U. We see that 
the discontinuity always decreases with increasing Av. In 
fact, the larger U is, the more abruptly it vanishes (on a 
scale of U) when Av > U. In this sense, the greater the 
asymmetry, the less discontinuous the energy derivative 
is, and the KS gap will be closer to the true gap. 

The situation is reversed when A = 1, as shown in Fig. 
31. Now the discontinuity grows with increasing Av. In 
this case, a large asymmetry puts the electron mostly on 
one site. When an infinitesimal of an electron is added, 
it goes to the same site, but paying an energy cost of U. 
On the other hand, if Av is small, the first electron is 
spread over both sites, and so is the added infinitesimal. 


































26 



Av/U 


FIG. 30 Derivative discontinuity as a function of Av for = 1, 
and = 5. 



Av/U 


FIG. 31 Derivative discontinuity for N — 1 as a, function of 
Av for = 1, and = 5. 


reducing the energy cost by a factor of 2. So Axe Ul2 
in the weakly correlated near-symmetric limit. 


C. Discontinuity around ni = 1 for iV = 2 

The derivative discontinuity manifests itself in many 
different aspects of DFT. We have already seen how it 
affects both energies and potentials as N is continuously 
moved across an integer. Here we explore how it appears 
even at fixed particle number, as correlations become 
strong. 

For our Hubbard dimer, with any finite Au, if U ^ Av, 
we know each rii is close to 1. The overwhelmingly large 
U localizes each electron on opposite sites. In the limit 
as U ^ oo, all fluctuations are suppressed, and the dimer 
becomes two separate systems of one electron each. For 
large but finite U, and finite Av, one is on the integer 
deficient side, and the other has slightly more than one 
electron. All the statements made above about Af passing 


through 2 now apply as rii passes through 1. 

We can see the effects in many of our earlier figures. In 
Fig. 7, the slope of F for f/ = 10 appears discontinuous 
at ni = 1. F contains the discontinuity in both Tg and 
Fxc in the limit U ^ oo. However, in reality, this curve 
is not really discontinuous. Zooming in on F near ni = 1, 
one sees that on a scale of 0{1/U), F is rounded. 


2 
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0 100 200 
Av 

FIG. 32 Plots of An and Am in HF as a function of Av 
for U = 100 (2t = 1), and the BALDA charge density. The 
crossover from the charge-transfer to the Mott-Hubbard regime 
happens at about U ~ Av. 

The classic manifestation already appears in Fig. 4, 
the occupation difference as a function of Av. To em¬ 
phasize the point, in Fig. 32, we plot several curves for 
U = 100. This is the discontinuous change from having 1 
particle on each site to 2 on one site that occurs. This is 
important because the common approximate density func¬ 
tionals miss this discontinuity effect. Explicit continuous 
functionals of the density cannot behave this way. For 
the SOFT case, this is embodied in the HF curves of Fig. 
9: No matter how strong the value of U, these curves are 
linear. In RHF, An versus Av never evolves the sudden 
step discussed above, as shown in Fig. 24. On the other 
hand, the BALDA approximation contains an explicit 
discontinuity at ni = 1 in its formulas, and so captures 
this effect, at least to leading-order in U. In this sense, 
both BALDA and UHF capture the most important effect 
of strong correlation. On the other hand, as discussed 
in Sec 7.C, UHF ‘cheats’, while BALDA retains the cor¬ 
rect spin singlet. If BALDA’s effects could be (legally) 
built into real-space approximations, they would be able 
to accurately dissociate molecules, overcoming perhaps 
approximate DFT’s greatest practical failure. 

However, in Fig. 33, we simply zoom in on the region 
of the plot near Av = U. In fact, the exact curve is S- 
shaped, with a finite curvature on the scale of t. Now we 
see that, although both UHF and BALDA reproduce the 
discontinuous effect, the details are not quite right. UHF 




/ 

/ 

/ 

/ 

/ 

. / . 


/ 

/ 

/ 

/ 

/ 

/ 

. / . 1 

/ 

/ : 
/ 

/ 

/ 

/ 




—Exact 
—RHF 
—UHF 
— BALDA 

. / . 

/ 

/ 

/ 

/ 

/ 


i 














































27 



95 100 105 

Av 


FIG. 33 Plots of An and Am in HF as a function of Av 
for U = 100 (2t = 1), and the BALDA charge density. The 
crossover from the charge-transfer to the Mott-Hubbard regime 
happens at about U ^ Av. 

is admirably close in shape to the accurate curve, but its 
slope is too great at ni = 1. BALDA is accurate to leading 
order in 1/t/, and captures beautifully the region Av a 
little larger than [/, but is quite inaccurate below that. 
The presence of the gap in the BALDA potentials leads to 
the incorrect discontinuous behavior near Av = 98. But 
once again we emphasize that the important feature is 
that these approximations do capture the dominant effect, 
and that BALDA does so without breaking symmetry. 


9. CONCLUSIONS AND DISCUSSION 

So, what can we learn from this exercise in applying 
DFT methods to the simplest strongly correlated system? 
Perhaps the most important point is that there is a large 
cultural difference between many-body approaches and 
DFT methodology, and a considerable barrier to commu¬ 
nication. In Sec. 3.C, we saw that even the definition 
of exchange is different in the two communities. The 
greatest misunderstandings come not from using different 
words for the same thing, but rather from using the same 
word for two different things. 

We can also see that the limitations of DFT calcula¬ 
tions are often misunderstood in the broader community. 
For example, the exact ground-state XC functional has a 
HOMO-LUMO gap that does not, in general, match the 
fundamental gap. The KS eigenvalues are not quasipar¬ 
ticle eigenvalues in general, and are in fact, much closer 
to optical excitations(A. Savin, 1998). Even the purpose 
of a DFT calculation is quite foreign to most solid-state 
physics. The modern art of DFT is aimed at producing 
extremely accurate (by physics standards) ground-state 
energies, and the many properties that can be extracted 


from those, rather than the response properties that are 
probed in most solid-state experiments, such as photoemis¬ 
sion. (Flipping the coin, most quantum chemists would 
never describe DFT energies as extremely accurate, as 
traditional quantum chemical ab initio methods are hyper 
accurate on this scale.) 

We also mention many aspects that we have not cov¬ 
ered here. For example, time-dependent DFT is based 
on a distinct theorem (the Runge-Gross theorem (Runge 
and Gross, 1984)), and provides approximate optical ex¬ 
citations for molecular systems (Burke et a/., 2005). The 
Mermin theorem(Mermin, 1965) generalizes the HK the¬ 
orem to thermal ensembles (Pribram-Jones et a/., 2014). 
There are many interesting features related to spin po¬ 
larization and dynamics, but very little is relevant to the 
system discussed here. There are also many non-DFT 
approaches, such as GW, which could be tested on the 
asymmetric dimer. 

We also take a moment to discuss how SOFT calcula¬ 
tions can be related to real-space DFT. One can easily add 
more orbitals to each site and create an extended Hub¬ 
bard model. For the H 2 molecule, adding just Pz orbitals 
and allowing them to scale yields a very accurate binding 
curve. But such an extension (beyond one basis function 
per site) is extremely problematic for SOFT(Harriman, 
1986), because it is no longer clear how to represent the 
‘density’. With 2 basis functions, should one use just the 
diagonal occupations, or include off-diagonal elements? In 
fact, neither one is satisfactory, as neither approaches the 
real-space density functional in the infinite basis limit. An 
underlying important point of DFT is that it is applied 
only to potentials that are diagonal in r, i.e., 'r’(r), and not 
diagonal in an arbitrary basis. This is a key requirement 
of the HK theorem, and is the reason why the one-body 
density n(r) is the corresponding variable on which to 
build the theory, and why the local density approximation 
is the starting point of all DFT approximations. 

This inability to go from SOFT calculations to real- 
space DFT calculations should be regarded as a major 
caveat for those using SOFT to explore DFT. Here we 
have shown many similarities in the behavior of SOFT 
functionals compared to real-space functionals. We have 
also proven some of the same basic theorems as those used 
in real-space DFT. But any results (especially unusual 
ones) that are found in SOFT calculations might not 
generalize to real-space DFT. The only way to be sure is 
to find a proof or calculation in real-space. On the other 
hand, SOFT calculations can be safely used to illustrate 
the basic physics behind real-space results (Stoudenmire 
et a/., 2012). 

Another limitation of SOFT can be seen already in 
our asymmetric Hubbard dimer. In a real heterogeneous 
diatomic molecule, say LiH with a pseudopotential for the 
core Li electrons, the values of U would be different on the 
two sites. But the basic DFT machinery only applies if 
the interaction is the same among all particles. And even 
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if it applies when both U and t become site-dependent, 
i.e., a one-to-one correspondence can be proven, it is 
unlikely that such studies would yield behavior that is 
even qualitatively similar to real-space DFT. 

Finally, we wish to emphasize once again the impor¬ 
tance of testing ideas on the asymmetric Hubbard dimer. 
Much (but not all) of the SOFT literature tests ideas on 
homogeneous cases. The essence of DFT is the creation 
of a universal functional, i.e., F[n] is the same no matter 
which specific problem you are trying to solve. The sym¬ 
metric case is very special in several ways, and there are 
no difficulties in applying any method to the asymmetric 
case. We hope that some of the results presented here 
will make that easier. 
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Appendix A: Exact solution, components, and limits 

In all appendices, we use dimensionless variables for 
brevity. Hence e = E/2t, u = I//2t, and u = Au/2t. 
Then, the energy of the singlet-ground-state is 

e= ^ (w-wsin(6»+^)) (Al) 

where 

w = \/3 [1 + z/^] + (A2) 

and 

cos(36>) = (9(z^^ — 1/2) + v?)u/w^. (A3) 

The coefficients of the minimizing wavefunction, Eq. 
(98), are 

= c , I3i^2 = c {u- e±u) , (A4) 

= 2 (z/^ + (e - uf (l -h e“^)) . (A5) 

The ground-state expectation values of the density 
difference and of the different pieces of the Hamiltonian 


are 


An = 4 z^ (n — e) 

(A6) 

V = —Au Anl2, 

(A7) 

II 

to 

1 

to 

T 

(A8) 

Fee = 4 y i U ((e - + y) . 

(A9) 


For fixed asymmetry z/, we can expand e in the weakly 
and strongly correlated limits: 

(i - (i++(1+ 

(AlO) 

where u = u/(l z^^)^/^. In the strongly correlated limit: 

e"' = -u-^ + (1 - -h 0{u-^). (All) 

We can also expand for fixed u around the symmetric 
limit: 

= -(u-r)+ (A12) 

2^ ^r{u + r) ^ ' 

where r = y/u‘^ 4. And the asymmetric limit: 

^asy ^ _^_^^_(2zy)-i_^^/2z/-2 + (i_4rt2^(2z/)-^ (Al3) 

Appendix B: Many limits of F(An) 

In this appendix we derive the limits that our param¬ 
eterization in Section 6 satisfies. Minimizing F of Eq. 
(104) with respect to we obtain a sextic equation for g: 

(4 + u^) + (y (3 + u^) - 1) y + 

2 y y+y (y (s + - (2 + y - 

2 «y(i-y)fl-y(i-y) = o (bi) 

where we define p = An/2. The solution defines gm{p)^ 
and F{p) = F{gm{p)^ p). Next we expand in several limits, 
and F\U^p] = F[I/, p, However, equation (Bl) can 
not be solved analytically in general. 


1. Expansions for g{p^u) 

We expand p in 4 different limits, which are built into 
po of Eq. (106) in Section 6. 

The weakly correlated limit corresponds to n ^ 1. We 
thus expand p(p, n) in powers of u for fixed p, 

CX) 

9{P, y = E u^/n\, (B2) 

n=0 

and insert the expansion into Eq. (Bl). The coefficients 
g^'^^ are found by canceling each term order by order in 
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Eq. (Bl), yielding 

= v/1 - P\ 9^^^ = 0, (B3) 

</'*’ = ^ (1 -( 1 + Ti>‘ - ^f'y 

Notice that ni ^2 = 1 i P so that to first order in [/, Eq. 
(104) yields the non-interacting kinetic energy functional 
ofEq. (43). 

For strongly correlated systems, we expand g in powers 
oil/u while holding p fixed 

CX3 

9{p,u) = '^ g^'"\p)u~^/n\, (B4) 

n=0 

and substitute back into Eq. (Bl) to find the coefficients. 
The result is 

5^°^ = y2/)(l -/>), = ly^, (B5) 

~(2) ^ ^i^p-^) ~(Q) 

8p 

Notice that this expansion breaks down at the symmetric 
point p = 0. 

The other kind of limit keeps u fixed. The symmetric 
limit is equivalent to p ^ 0. We expand g in powers of p 
while holding u fixed. 

CX3 

9 {p,u) = /n\, (B6) 

n=0 


compute Ec in each regime: 

ec = -5 + uh{g, p) - | (1 + p) + \/l - p. 

where h{g^p) is defined in Eq. (105). Then, as rt ^ 0, 
Cc ^ where 

e:?(p = -y(i-P)®/^(i-«P\/r^). (BIO) 

Similarly, as ^ oc, Cc ^ ^ where 

= ~2 + “ P (\/l + P ~ V^)— 

(Bll) 

An alternative expansion is to fix u and expand in p. As 
p ^ 0, Cc ^ where 

e^y^^ip) = 1 - ^l+d)' (B12) 



As p ^ 1, Cc ^ where 

p/2 (^-^+u^y (B13) 

where p = 1 — p. 

3. Order of limits 


and substitute back into Eq. (Bl) to find the coefficients. 
The result is 


(B7) 

where r = y/l + {u/2)‘^. 

The asymmetric limit is equivalent to p ^ 1. We 
expand g in powers of p = 1 — p for fixed u: 

oo 

n=0 




5 ( 2 ) = 

^ 2 


0 u^l2{u^l2 

H- -——— 


and substitute back into Eq. (Bl). The result is 

pi/2) = yp2, p3/2) = -3pi/2)/8 (B9) 

= 12 



2. Limits of the correlation energy functional 

Now that we have expressions for g in all four limits 
we can use our expression for F, eq. (104), Tg, and t/n to 


Finally, we look at how these expressions behave when 
both parameters are extreme. The weakly correlated limit 
has no difficulties near the symmetric point: 


e-(p ^ 0) = e%y^{u ^ 0) 



V? p^ 


(B14) 


In the asymmetric limit, there are also no problems: 


e-(p ^ 1) = ^ 0) 

= M2p/2 , (B15) 

Thus, the expansion in powers of u is well-behaved, and 
there are no difficulties using it for sufficiently small u. 
In the symmetric case, one sees explicitly that the radius 
of convergence of the expansion is u = 2. 

On the other hand, the strong coupling limit is more 
problematic. Expanding the strong-couping functional 
around the symmetric limit, we find 


(P ^ 0) 


u 1 

■2+^“4^ 


\/^ + P 



Au) ’ 
(B16) 
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while reversing the order of limits yields: 


psym 


oo) = + 1 - ^- ^ ~ ^ ~ ^ 


2 

(B17) 

Note the difference beginning in the third terms, i.e., at 
first-order in 1/rt, even for p = 0. Thus for the Hubbard 
dimer, approximations based on expansions around the 
strong-coupling limit are likely to fail for some values of 
the density. 


Appendix C: Proofs of Energy Relations 

Using the notation established in Section 6, we prove 
some simple relations about the energy and its compo¬ 
nents. Start with the general expression for the energy, 
Eq. (104) and (105), 

e = min [-g + uh{g, p) - vp]. (Cl) 

P,9 

First take p —> 0. The second term reduces to 
u (^1 — ^/l — /2. Then let g ^ 0, resulting in h ^ 0. 

This yields e ^ 0 and therefore the exact e < 0. This 
process corresponds to choosing a trial wavefunction, and 
by Rayleigh-Ritz, the ground-state wavefunction will pro¬ 
duce a value equal to or below the trial result. 

In Hartree-Fock, g reduces to ^hf = Then, 

= mme(gHF(p),p) > e- (C2) 

P 


The first term is less than zero by definition but the second 
needs more unraveling. To begin, from Eq. (104), 


dl 

dg 


dh 

dg 


so, at the solution 


dh 1 
dg u 


(C7) 


(C8) 


For A near 1, Suppose g{X) g{l) + (A — l)^'(l), and 
expand dh/dg\g(^x) in g{X) around ^(1): 



dh 

dg 


+ (A-1)^'(1) 

^(1) 


d^h 

dg‘^ 


9d) 


(C9) 


The first term on the left is l/{Xu) ^ {2 — X)/u. After 
some algebra. 


dg 


dX 


A=1 



(CIO) 


Since the hopping term of / is linear in d^f/dg^ = 
d^h/dg^. The energy is a minimum at g so d^f/dg^ > 0, 
thus dg/dX > 0. Together, this results in 


dU^/dA<0, (Cll) 


the adiabatic connection integrand is monotonically de¬ 
creasing as seen in Fig. 21. 


This shows that = e — > 0, as in Fig. (5). The 

minimization can be performed analytically though it 
involves solving the quartic polynomial 


— , ^ u p — u = {). 


(C3) 


Similarly, a DFT exact exchange (EXX) calculation is 
defined by 


Appendix D: BALDA Derivation 

For an infinite homogeneous Hubbard chain of density 
n = 1 + X, the energy per site (in units of 2 1) is given 
approximately by 

= ux0{x)a{x^ P(U)) (Dl) 


= e(5HF(/3m), Pm) > min e(5HF(p), p) (C4) 

P 

where is the minimizing density for the many-body 
problem. This yields = e — and > 

(Gross et a/., 1996). 

For the kinetic energy alone, t = —g{pm)^ and 


where 0{x) is the Heaviside function and 

a{x,f3) = —— sin(7r(l — \x\)/f3) /tt. (D2) 

TT 

The function /3{u) varies smoothly from 1 at = 0 to 2 
diS u ^ oo(Lima et a/., 2003), and satisfies 


ts = min [-g{p)] = -^/l- p^. (C5) 


This results in tc > 0 since the KS occupation difference 
is defined to minimize the hopping energy. This combined 
with the above implies Uc < 0, as in Eq. (79). 

For the adiabatic connection integrand, take a deriva¬ 
tive of Eq. (Ill): 


du^ ^ Uc{p,X) . dhdg 
dX ~ A dgdX' 


(C6) 


Q!(0, /3) = —4 



MO MO 

^ [1 + exp(u^))] 


(D3) 


This simple result is exact as rt ^ 0, rt ^ oc and at 
n = 1, and a good approximation (accurate to within a few 
percent) elsewhere (Lima et a/., 2003) to the exact solution 
via Bethe ansatz(Lieb and Wu, 1968). In principle, /3 
depends on n, and this dependence has been fit in later 
work(Franca et al.^ 2012b). Here, we use the simpler 
original version of a function of u only. In fact, the 
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solution to Eq. (D3) can be accurately fit (error below 
1%) with a simple rational function, 




2 au^ hv? 
1 + C14 + hv? 


(D4) 


with coefficients a = 2c —7r/4 and b = {a —c)/ log 2 chosen 
to recover the small-i4 behavior to first-order, and the large 
u behavior to first order in 1/u^ and c = 1.197963 is fit to 
P{u). This is useful for quick implementation of BALDA. 

At 14 = 0, the hopping energy per site is just 

=-sin(7r(l - |a:|))/7r, (D5) 


while the Hartree-exchange energy per site is a simple 
local function: 


<xV = un^/A. 


(D6) 


Thus the correlation energy per site is just 


~unif ~unif pinif ^~unif 

to c 6g . 


(D7) 


The BALDA approximation is then 

eBALDA ^ ~u„if ^ ^ ~u„if ( 08 ) 


Since the exchange is local, BALDA is exact for that 
contribution, and only correlation is approximated. Since 
ni ^2 = 1± An/2, x = ±An/2 for sites 1 and 2 respectively. 
The BALDA HXC energy is then: 

Crxc= -2 (a(An/2, U) - 0)) + w| An|/2, 

(D9) 

and was inserted into the KS equations (Sec 3.C) to find 
the results of Sec 7.B. 


Appendix E: Mean-Field Derivation 


The ME hamiltonian for the Hubbard dimer can be 
written in the number basis |lcr, 2g) as follows 


H, 


MF 


-/^vf -t \ 
-t Avf ) 


(El) 


with cr = ±1 for spin up and down respectively. Setting 
M = mi +7742 and A" = ni +n 2 as the total magnetization 
and particle number of the system, the eigenvalues are 

TT /eff 

(E2) 

tf = 2t^{Av%«/2tY + I, 

Av®® = Av — ^ (An — (j Am). 

The total energy of the system is 

=e_,t + e+,t-[/H (E3) 

EAFm = e_,t + e_4 - [/«, (E4) 


where the Hartree term is written as 

TT U ^ 

= - (AT^ - + An^ - Am^). (E5) 

8 


Depending on whether is larger or smaller than 

the ground-state of the system may be ferromag¬ 
netic (N = 2, \M\ = 2) or antiferromagnetic (N = 2, 
M = 0, |A777| > 0). The paramagnetic state is a specific 
case of the AEM state with Am = 0. Explicitly, for the 
ferromagnetic state we have the eigenstate energies and 
self-consistency equations 


An = Am = Av/ \/4+ Av‘^ (E6) 

ezF,t = tVaFF^/2 (E7) 

On the other hand, the M = 0 state (| AttiI > 0 is AEM, 
Am = 0 is PM) corresponds to the eigenvalues, 

e-,t = {U- tf)/2, e_4 = {U - tf)/2, (E8) 

and self-consistency equations 


An = E 


Av^^ 


Am = a 




(E9) 


and the expressions for Av^^ and are given in Eq. 
(E2). The self-consistency procedure needs to be carried 
out numerically in this case. 

The total energy can also be written as 


^AFM,PM ^ ^ 41 _ A77^ - Am‘^ 

2 ^ 4 ^ 


tfEtf 


(ElO) 


In the PM case, the expressions can be simplified to give 

2 Av — U An 


An = 


ij {Av — U An/2)^ + 4 
for the occupations and 


(Ell) 



(E12) 


Appendix F: Relation between Hubbard model and 
Real-space 


To show how SOFT and real-space DFT are connected, 
begin with the one-electron dimer, Hj, with the protons 
separated by R. Use a basis of the exact atomic Is orbitals, 
one on each site. This is a minimal basis in quantum 
chemistry. Then 



1 

|r — i7z| 


(El) 
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where the bond is along the 2 ;-axis. Then the matrix 
elements of h in the basis set of atomic orbitals are: 

vi = V2 = eA j{R)^ t = s{R) €a + k{R) (F2) 
where is the atomic energy (- one Rydberg here) and 
s{R) = {A\B) = e-^{l + R + R^/3) 
j{R) = = -(i/i? - + l/R)) 

= = + (F3) 

yielding the textbook eigenvalues (for the generalized 
eigenvalue problem): 

+ {j ± k)/ (1 ± s). (F4) 

Of course, the orbitals can always be symmetrically or- 
thogonalized in advance (Lowdin, 1950), in which case 

I'ortho = eA + {-j + ks)/ (s^ - 1), (F5) 

iortho = -(sj - k)/(s^ - 1). (F6) 

Although physics textbooks often set the overlap to zero, 
this is inconsistent, as the size of the overlap is comparable 
to k(R), say. Setting the on-site potential to zero (but re¬ 
adding its value to the energy) and using tortho? makes the 
solution Eq. (29) of the text produce the exact electronic 
energy in this minimal basis. 

But quantum chemistry textbooks note that this calcu¬ 
lation is horribly inaccurate, yielding a bond-length of 2.5 
Bohr and a well depth of 2.75 eV. Inclusion of a orbital 
on each site, and allowing the lengthscale of each orbital 
to vary, produces almost exact results of 2.00 Bohr and 
4.76 eV. Thus, even in this simple case, more than one 
orbital per site is needed to converge to the real-space 
limit. 

Next we consider repeating the minimal-basis calcula¬ 
tion with one nuclear charge replaced by value Z. This 
yields an asymmetric tight-binding problem for which 
the orbitals can be orthogonalized and values of Av and 
t deduced as a function of R. But note that changing 
Z will change both Av and t simultaneously, unlike our 
asymmetric SOFT dimer, where only Av changes. In 
real-space DFT, the kinetic energy functional remains the 
same, of Eq. (20), for all R and every Z. 

The situation is even more complicated for H 2 and its 
asymmetric variants. Clearly U becomes a function of 
i7, but there are also several independent off-diagonal 
matrix elements that are R dependent. Again, all change 
as a function of both R and Z, but none of this occurs in 
SOET. In real-space DET, Tg is still the von Weisacker 
functional, Uh always the Hartree energy, and the exact 
^xc['^] is independent of R and Z, but always produces 
the exact energy when iterated in the KS equations. 
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